Working Wonders with ADPL Math - Ep 03: Data Reduction Fundamentals - Page 2 of 3

Compression using SVD: a minimal working example

Say we have a nx3 matrix «M», for which the first two columns are linearly independent (geometrically speaking, those vectors are «non parrallel»), and the third one is a linear combination of the first two.

From linear algebra, we know that this matrix has rank 2, and hence is non-invertible. So far so good. But that’s not our point, since we don’t want to solve for the matrix. Quite the contrary, we assume that there are redundancies (or repetitions) in its vectors, so that it’s non invertible, hence amenable to compression. Then, from this perspective, we know that it would be more economical to describe this matrix as «its first vector is like this, the second vector like that, and the third is a linear combination of the first two, with the following coefficients».

Well, SVD does exactly that. Let’s see how.

Mathematically, SVD will factor the original matrix M (nxm) into a product of 3 matrices, namely

$M=U \Sigma V$

U=left singular vector matrix (nxn)
$\Sigma$ = diagonal matrix, containing singular values (nxm)
V’=right singular vector matrix (mxm)

In this factorization, the number of non-zero entries in is equal to the rank of M. Obviously, if M is rectangular, M has rank at most equal to min(n,m). But here, we are specifically looking at situations where redundancies exist, so we expect the rank r to be much lower than that.

Right, let’s see how it works within ANSYS. Assume for simplicity a 3×3 matrix.

First column vector $M_1= \begin{bmatrix} 1 \\1 \\ 1 \end{bmatrix}$

Second column vector $M_2= \begin{bmatrix} 1 \\ -1 \\ 1 \end{bmatrix}$

Third column vector $M_3=M_1+M_2$

The corresponding input commands would be:

n=3	! nb of "DOFs"
m=3	! nb of "snapshots"

! Generate dummy "snapshots" matrix
*DMAT,dummyMat,D,ALLOC,n,m

! First and second vectors (should be linearly independent, hence rank=2)
! first column 
dummyMat(1,1)=1	
dummyMat(2,1)=1	 
dummyMat(3,1)=1	 
! second column
dummyMat(1,2)=1		
dummyMat(2,2)=-1
dummyMat(3,2)=1
! Add a third column, a linear combination of the first two
alpha=1
beta=1
dummyMat(1,3)=alpha*dummyMat(1,1)+beta*dummyMat(1,2)
dummyMat(2,3)=alpha*dummyMat(2,1)+beta*dummyMat(2,2)
dummyMat(3,3)=alpha*dummyMat(3,1)+beta*dummyMat(3,2)

As expected, we get the following

DUMMYMAT: 
 [1,1]: 1.000e+00 [1,2]: 1.000e+00 [1,3]: 2.000e+00 
 [2,1]: 1.000e+00 [2,2]:-1.000e+00 [2,3]: 0.000e+00 
 [3,1]: 1.000e+00 [3,2]: 1.000e+00 [3,3]: 2.000e+00

From there, we can simply invoke the matrix compression by first allocating for the Matrix U, which will contain the left singular vectors, corresponding to singular values above a user defined threshold (more on this later). To do that, we simply copy the original matrix (it will be read and replaced by a skinnier version, same number of rows, but lower number of columns). The sigma matrix is diagonal and will be stored in a vector. Note that both vector sigmaVec and matrix Vtran will be created on the fly, and don’t require manual allocation.

! Perform Singular Value Decomposition 
*DMAT,U,D,COPY,dummyMat	

! M=U SIGMA V*
threshold=1e-2
*COMPRESS,U,SVD,threshold,SigmaVec,Vconj

The output reads as follows:

 *COM COMMAND : USING SVD ALGORITHM  Allocate a [3] Vector : SIGMAVEC

  Allocate a [3][3] Dense Matrix : VCONJ
 (TOLERANCE = 0.01)
                MATRIX SIZE AFTER COMPRESSION : 2

 U: 
 [1,1]:-7.071e-01 [1,2]: 1.504e-16 
 [2,1]: 7.673e-17 [2,2]:-1.000e+00 
 [3,1]:-7.071e-01 [3,2]: 5.395e-17 
 SIGMAVEC :
 Size : 3
  3.464e+00   1.414e+00   7.009e-17  

 VCONJ: 
 [1,1]:-4.082e-01 [1,2]:-4.082e-01 [1,3]:-8.165e-01 
 [2,1]:-7.071e-01 [2,2]: 7.071e-01 [2,3]:-7.477e-17 
 [3,1]:-5.774e-01 [3,2]:-5.774e-01 [3,3]: 5.774e-01

A number of important features are already visible:

– the U matrix has been compressed to size 3×2, consistent with the fact that only 2 singular values are above the threshold (the last one being a numerical zero)

– U and V are unitary matrices (vectors have unit norm, and are orthogonal), hence U does not contain our «basis» vectors, but rather an orthogonal basis for the r-dimensional space in which the matrix M projects the n-dimensional space.

– Vconj has 3 columns, but also 3 rows: however, the last row corresponds to the amplitudes of the discarded 3^rd right singular vector, and has no usage.

This decomposition and the «padding» is explained in an article by Gregory Undersen, see [2]. In order to check that everything ran as expected, we might want to re-generate the original matrix from its truncated factors, as follows:

! Number of Singular values above threshold
nbSingularValues=U_coldim	

! CHECK: REGENERATE ORIGINAL MATRIX FROM COMPRESSED 
! Step1: Generate SigmaMat (diagonal matrix) from vector SigmaVec
*SMAT,SigmaMat,D,ALLOC,DIAG,nbSingularValues
*do,ind,1,nbSingularValues
	SigmaMat(ind,ind)=SigmaVec(ind)
*enddo

! STEP2: GENERATE TRUNCATED RIGHT VECTORS MATRICE 
*DMAT,resizedVconj,D,COPY,Vconj
*DMAT,resizedVconj,D,RESIZE,nbSingularValues,nbSnapshots

!*MULT, M1, T1, M2, T2, M3 (M3=M1(T1)*M2(T2))
*MULT,SigmaMat,,Vconj,,SigmaVconj
*PRINT,SigmaVconj
*MULT,U,,SigmaVconj,,shouldBeDummyMat
*PRINT,shouldBeDummyMat

The output is consistent with our expectations:

SHOULDBEDUMMYMAT: 
 [1,1]: 1.000e+00 [1,2]: 1.000e+00 [1,3]: 2.000e+00 
 [2,1]: 1.000e+00 [2,2]:-1.000e+00 [2,3]:-1.113e-16 
 [3,1]: 1.000e+00 [3,2]: 1.000e+00 [3,3]: 2.000e+00

This is all well, but one could (and in fact, should) argue that this particular example was mathematically designed to showcase repeated vectors identification, and obviously works to perfection. In a real situation, one should expect a gradual decrease of singular values or, stated differently, a fuzzier boundary between «needed» and «superfluous» components, which renders the singular values selection much more empiric.

A bulletproof -though naïve- approach would be to iteratively generate data using an increasingly large number of singular vectors, until a satisfactory result is obtained, but this may not even be required. As it appears, since U and V matrices are unitary, so that only the singular values contain information about amplitude (i.e. they are scaling factors). What’s more, those values are provided in decreasing order, and the total energy of the original signals (the square of the Frobenius norm of the original matrix) is equal to the sum of the squared singular values.

Luckily for us ANSYS dutifully follows this rule. By construction, the Froebenius (energy) norm of our dummy matrix M is the sum of the norm of its vectors, namely:

$||M_1||^2=1^2+1^2+1^2=3$

$||M_2||^2=1^2+(-1)^2+1^2=3$

and

$||M_3||^2=2^2+(0)^2+2^2=8$

So that its total energy norm is :

$||M||=3+3+8=14$

This value compares favorably with the sum of the squared singular values calculated by ANSYS:

$Tr(\Sigma \Sigma ^t)=\sigma_1^2+\sigma_2^2=(3.464)^2+(1.414)^2=14$

Hence, by selecting a number of singular values, the amount of the total energy is known beforehand, and the energy loss (sum of squared unused singular values) can be controlled. To reemphasize, compression will only be efficient when the singular values decrease sufficiently rapidly, which may seem difficult to detect with the naked eye. A very clear explanation has been given by Pr.Alex Townsend in [3], which sheds light on the fact that propagating waves phenomena are ill-suited for SVD compression.