Eigen Pairs

In this section, let us develop our understanding of linear transformation by breaking it down into simpler linear transformations. 

  • Let T be a linear transformation from \mathbb {R}^n to \mathbb {R}^n.

Suppose:

  • B is a basis of \mathbb {R}^n
  • V is the span of some of the vectors in B
  • W is the span of the remaining vectors in B.

Then:

  • Any vector in \mathbb {R}^n can be written as the sum of a vector v in V and a vector w in W.
  • Since T(v + w) = T(v) + T(w), we can see how T behaves on all of \mathbb {R}^n if we know how it behaves on V and on W.
  • This decomposition is particularly helpful if V and W are chosen so that T behaves in a simple way on V and on W.

Given such a decomposition of \mathbb {R}^n into the vector spaces V and W, the same idea can be applied to split V and W into lower-dimensional vector spaces and repeat until no more splits are possible. The most optimistic outcome of this procedure would be that we get all the way down to n one-dimensional subspaces and that T acts on each of these subspaces by simply scaling each vector in that subspace by some factor. In other words, we would like to find n vectors v for which T(v) is a scalar multiple of v. From this, the following definition can be obtained. 

Mathematical Definition:

  • An eigenvector v of an n \times n matrix A has the property Av = \lambda v for some \lambda \mathbb {R}.
  • \lambda is the eigenvalue of A.
  • The eigenvector together with its eigenvalue is called eigenpair.

Being more general, an eigenvector is a vector whose direction remains unchanged when a linear transformation is applied to it.

In the figure above, the red vectors are the eigenvectors because their directions didn’t change even after the application of a linear transformation(scaling in this case) while the magenta vector is not an eigenvector as its direction is changed after the linear transformation.

Q: Find the eigen pairs of a 3 \times 3 matrix A = \begin{bmatrix} 1.5 & 0 & 1\\-0.5 & 0.5 & -0.5 \\-0.5 & 0 & 0 \end{bmatrix}

				
					%use s2

// Create a matrix
val A = DenseMatrix(arrayOf
        (doubleArrayOf(1.5, 0.0, 1.0), 
         doubleArrayOf(-0.5, 0.5, -0.5), 
         doubleArrayOf(-0.5, 0.0, 0.0)))
// eigen decomposition
val eigen = Eigen(A, Eigen.Method.QR, 0.0)
println("eigenpairs:")
for (i in 0 until eigen.size()) {
    // get the i-th eigenvalue
    val eigenvalue = eigen.getEigenvalue(i)
    // get the property associated with the eigenvalue
    val property = eigen.getProperty(eigenvalue)
    // get the eigen vector from property
    val eigenVector = property.eigenVector()
    println("eigenvalue: $eigenvalue; eigenvector: $eigenVector")
}
				
			
eigenpairs:
eigenvalue: 1.0; eigenvector: [-2.000000, 1.000000, 1.000000] 
eigenvalue: 0.5; eigenvector: [-0.000000, 1.000000, 0.000000] 
eigenvalue: 0.49999999999999994; eigenvector: [-0.000000, 1.000000, 0.000000] 

Note: Every non-zero vector is an eigenvector of the identity matrix with eigenvalue 1.

Eigen Decomposition

Let’s say, A is an n \times n having linearly y independent eigenvectors, then one-dimensional subspaces spanned by each of these vectors can be thought of as (not necessarily orthogonal) axes along which A acts by scaling.

Speaking in matrix terms, we can define Q to be the matrix with the eigenvectors of A as columns. Then from the definition of an eigenpair, we have:

AQ = QD

where  D is a diagonal matrix whose diagonal entries are the eigenvalues (in the order corresponding to the columns of Q) and whose other entries are zero.

Thus, we can conclude that A = QDQ^{-1}, where  D is a diagonal matrix and A is diagonalizable.

S2 makes it easier to find the diagonal matrix D for a given diagonalizable matrix A.

				
					%use s2

// Create a symmetric matrix
val A = DenseMatrix(SymmetricMatrix(
    arrayOf(doubleArrayOf(1.0),
            doubleArrayOf(2.0, 3.0), 
            doubleArrayOf(4.0, 5.0, 6.0), 
            doubleArrayOf(7.0, 8.0, 9.0, 10.0))))
// eigen decomposition
val precision = 1e-4
val decomp = EigenDecomposition(A, precision)
val D = decomp.D() // get the diagonal matrix which contains the eigen values
println("eigen decomposition result: $D")

// verification
val Q = decomp.Q();
val Qt = decomp.Qt();
val B = Q.multiply(D).multiply(Qt); // QDQ' = A = B
println("$A =\n$B")
				
			
eigen decomposition result: 4x4
	[,1] [,2] [,3] [,4] 
[1,] 24.638176, 0.000000, 0.000000, 0.000000, 
[2,] 0.000000, -0.038687, 0.000000, 0.000000, 
[3,] 0.000000, 0.000000, -0.513527, 0.000000, 
[4,] 0.000000, 0.000000, 0.000000, -4.085962, 
4x4
	[,1] [,2] [,3] [,4] 
[1,] 1.000000, 2.000000, 4.000000, 7.000000, 
[2,] 2.000000, 3.000000, 5.000000, 8.000000, 
[3,] 4.000000, 5.000000, 6.000000, 9.000000, 
[4,] 7.000000, 8.000000, 9.000000, 10.000000,  =
4x4
	[,1] [,2] [,3] [,4] 
[1,] 1.000000, 2.000000, 4.000000, 7.000000, 
[2,] 2.000000, 3.000000, 5.000000, 8.000000, 
[3,] 4.000000, 5.000000, 6.000000, 9.000000, 
[4,] 7.000000, 8.000000, 9.000000, 10.000000, 

Positive Definite Matrices

 A symmetric matrix A is called a positive definite matrix if its eigenvalues are all positive. A symmetric matric matrix A is a positive semidefinite matrix if its eigenvalues are all non-negative.

Equivalently, a matrix A is positive semidefinite if x'Ax \geq 0 for all x.

Analogously, a symmetric matrix A is called a negative definite matrix if its eigenvalues are all negative. A symmetric matric matrix A is a negative semidefinite matrix if its eigenvalues are all non-positive.

In S2, we have a constructor “PositiveDefiniteMatrixByPositiveDiagonal” that converts a matrix into a symmetric, positive definite matrix, if it is not already so, by forcing the diagonal entries in the eigen decomposition to a small non-negative number.

				
					%use s2
// define a matrix
var A = DenseMatrix(arrayOf(
    doubleArrayOf(1.0, 2.0, 3.0),
    doubleArrayOf(4.0, 5.0, 6.0),
    doubleArrayOf(-1.0, -2.0, -3.0)))

// A_PDM = Positive Definite Matrix of A
val epsilon = 1e-15
val A_PDM = PositiveDefiniteMatrixByPositiveDiagonal(A,epsilon,1.0)
println("POSITIVE DEFINITE MATRIX OF\n\n $A \n\n IS \n\n $A_PDM\n")
				
			
POSITIVE DEFINITE MATRIX OF

 3x3
	[,1] [,2] [,3] 
[1,] 1.000000, 2.000000, 3.000000, 
[2,] 4.000000, 5.000000, 6.000000, 
[3,] -1.000000, -2.000000, -3.000000,  

 IS 

 3x3
	[,1] [,2] [,3] 
[1,] 2.286186, 2.413588, 0.605269, 
[2,] 2.413588, 5.529211, 1.135817, 
[3,] 0.605269, 1.135817, 1.284835, 

We also have another constructor in S2 named “PositiveSemiDefiniteMatrixNonNegativeDiagonal” which converts a matrix into a symmetric, positive semi-definite matrix, if it is not already so, by forcing the negative diagonal entries in the eigen decomposition to 0.

				
					%use s2
// define a matrix
var A = DenseMatrix(arrayOf(
    doubleArrayOf(1.0, 2.0, 3.0),
    doubleArrayOf(4.0, 5.0, 6.0),
    doubleArrayOf(-1.0, -2.0, -3.0)))

// A_PSDM = Positive Semi-Definite Matrix of A
val epsilon = 100.0
val A_PSDM = PositiveSemiDefiniteMatrixNonNegativeDiagonal(A,epsilon)
println("POSITIVE SEMI-DEFINITE MATRIX OF\n\n $A \n\n IS \n\n $A_PSDM\n")
				
			
POSITIVE SEMI-DEFINITE MATRIX OF

 3x3
	[,1] [,2] [,3] 
[1,] 1.000000, 2.000000, 3.000000, 
[2,] 4.000000, 5.000000, 6.000000, 
[3,] -1.000000, -2.000000, -3.000000,  

 IS 

 3x3
	[,1] [,2] [,3] 
[1,] 7.100233, -0.000000, -0.000000, 
[2,] -0.000000, 7.100233, -0.000000, 
[3,] -0.000000, -0.000000, 7.100233, 

Spectral Theorem:

If A is an n \times n symmetric matrix, then A is orthogonally diagonalizable, which means A has n eigenvectors which are pairwise orthogonal.

Conversely, every orthogonally diagonalizable matrix is symmetric.

Polar Decomposition

Before proceeding with the definition of polar decomposition, let us make ourselves familiar with “Gram Matrix”:

Gram Matrix: If A is an m \times n matrix, then A'A is its Gram matrix. The Gram matrix of A is always positive semidefinite.

Q: Determine the Gram Matrix of A = \begin{bmatrix}1 & 2 & 3\\4 & 5 & 6\\-1 & -2 & -3\end{bmatrix}.  

				
					%use s2
// define a matrix
var A = DenseMatrix(arrayOf(
    doubleArrayOf(1.0, 2.0, 3.0),
    doubleArrayOf(4.0, 5.0, 6.0),
    doubleArrayOf(-1.0, -2.0, -3.0)))

// A_T = Tranpose of A
val A_T = A.t()

//A_T*A = Gram Matrix of A
val A_GM = A_T.multiply(A)
println("THE GRAM MATRIX OF \n $A \n\n IS \n\n $A_GM\n")
				
			
THE GRAM MATRIX OF 
 3x3
	[,1] [,2] [,3] 
[1,] 1.000000, 2.000000, 3.000000, 
[2,] 4.000000, 5.000000, 6.000000, 
[3,] -1.000000, -2.000000, -3.000000,  

 IS 

 3x3
	[,1] [,2] [,3] 
[1,] 18.000000, 24.000000, 30.000000, 
[2,] 24.000000, 33.000000, 42.000000, 
[3,] 30.000000, 42.000000, 54.000000, 

Let’s define matrix  \sqrt{A'A} to be QD^{1/2}Q', where QDQ' is the diagonalization of A'A and D^{1/2} is the matrix obtained by taking the diagonal entries of D. Then  \sqrt{A'A} is symmetric and satisfies:

 \sqrt{A'A} \sqrt{A'A} = QD^{1/2}Q' QD^{1/2}Q' = A'A

The matrix  \sqrt{A'A} is easier to understand than A because it is symmetric and positive definite. However, it transforms space in nearly same way as A. If  x \in \mathbb {R}^n, then:

 \left|Ax\right|^2 = x'A' Ax = x' \sqrt{A'A} \sqrt{A'A} x = \left| \sqrt{A'A}x \right|^2

In other words, for all x, the images of x under A and under  \sqrt{A'A} have equal norm. This means that for each x \in \mathbb {R}^n, there is an orthogonal transformation from the range of  \sqrt{A'A} to the range of A which sends
Ax to  \sqrt{A'A}x. That means, this orthogonal transformation is the same for all x.

Fig: From the above figure, we can observe that the grid-line images under A and  \sqrt{A'A} are of the same shape.

 \implies They are related by orthogonal transformation.

This orthogonal transformation can be extended to an orthogonal transformation on  \mathbb {R}^n even if the range of  \sqrt{A'A} is not all of  \mathbb{R}^n. Hence, we arrive at a polar decomposition.

Polar Decomposition:

For any n \times n matrix A, there exists an orthogonal matrix R such that

A = R \sqrt{A'A}

 

Let us implement the polar decomposition using S2.

				
					%use s2
//define matrix A
var A = DenseMatrix(arrayOf(
    doubleArrayOf(4.0, 5.0),
    doubleArrayOf(6.0, 7.0)))
 
//AT is transpose of A
var AT = A.t()
//ATA = AT*A
var ATA = AT.multiply(A)
println(ATA)
				
			
2x2
	[,1] [,2] 
[1,] 52.000000, 62.000000, 
[2,] 62.000000, 74.000000,
				
					%use s2
// define a matrix sqrt(A'A) from previous results
var RATA = DenseMatrix(arrayOf(
    doubleArrayOf((52.0).pow(0.5), (62.0).pow(0.5)),
    doubleArrayOf((62.0).pow(0.5), (74.0).pow(0.5))))

//Matrix A
var A = DenseMatrix(arrayOf(
    doubleArrayOf(4.0, 5.0),
    doubleArrayOf(6.0, 7.0)))

//calculate inverse of sqrt(A'A)
val INVRATA = Inverse(RATA)

// R = A*Inverse(sqrt(A'A))
val R = A.multiply(INVRATA)
println("R =\n\n $R\n")

//verification
var B = R.multiply(RATA)
println("\nR*sqrt(A'A) = \n\n $B \n\nSAME AS\n\n A = \n\n$A\n")
				
			
R =

 2x2
	[,1] [,2] 
[1,] -153.822883, 141.380679, 
[2,] -108.655461, 100.269860, 


R*sqrt(A'A) = 

 2x2
	[,1] [,2] 
[1,] 4.000000, 5.000000, 
[2,] 6.000000, 7.000000,  

SAME AS

 A = 

2x2
	[,1] [,2] 
[1,] 4.000000, 5.000000, 
[2,] 6.000000, 7.000000,