principal axis transformation

Principal Axis Transformation and SVD PhL 6.9.12 Motivation: this is a similar "active vs passive" topic to the vector discussion nearby, but the objects of interest here are matrices instead of vectors. I have never really gotten this subject written down. It occurs with the inertia tensor, the stress tensor, the Ng covariance tensor, and so on. Standard tensor notation and graphical representation. Use developmental tensor doc notation. Here then is how rank 1 and rank 2 tensors transform : V'a = RaαVα // contravariant vector V' = R V // previous line in vector notation x' = R x // example of a vector transformation M'ab = Raα RbβMαβ // contravariant rank 2 tensor M' = RMRT = RMR-1 // previous line in matrix notation So the above equations say that V and x are vectors with respect to transformation R, and that M is a rank-2 tensor with respect to R. Notice that the idea that V' = R V works with two of my three vector viewing pictures: Now consider another set of equations '(i) = R-1 (i) // previous line applied to ith basis vector u(i) '(i)a = R-1ab (i)b // previous line in components Here we see that (i) are vectors with respect to transformation R-1 (not R as for the vectors in the first list of equations above). If we multiply both sides of the last equation above by (j)a we get (j)a '(i)a = (j)a R-1ab (i)b = R-1ji = Rij But we also know that (j)a '(i)a = (j) '(i) and (j)a '(i)a = δj,a'(i)a = '(i)j To summarize, we have (j) '(i) = '(i)j = Rij The fact that '(i)j = Rij says that row i of matrix Rij is in fact the vector '(i), so one could write RT = ['(1), '(2), .... '(n)] . Once we have the above two sets of equations, we can explore things. One thing we find is that V = Vi (i) = V'i '(i) is consistent. Proof: Write the right equality as Vi (i)a = V'i '(i)a Then consider: V'a = RaαVα '(i)a = R-1aβ(i)β Then we find that V'i '(i)a = (RiαVα) (R-1aβ(i)β) = (RiαVα) (R-1aβδi,β) = (RiαVα) (R-1ai) = (RiαVα) (RTai) = (RTaiRiα)Vα = δa,α Vα = Va Vi (i)a = Vi δi,a = Va Now things are compatible with the passive view picture V = Vi (i) = V'i '(i) // names of unit vectors don't match picture but OK So when we "view" the vector V from the frame S', we find that it has certain components V'i which are different from the components Vi that we "view" from frame S. Focus on the rank 2 tensor. Now how do we do this "view" stuff for a rank 2 tensor? M'ab = Raα RbβMαβ Well, we can expand M in two different ways, following tensor doc M = Mij (i) (j)T analogous to V = Vi (i) M = M'ij '(i) '(j)T analogous to V = V'i '(i) I am sure that I could verify that both lines are consistent. Here we go Mab = Mij (i)a (j)Tb = Mij δi,aδj,b = Mab M'ij '(i)a '(j)Tb = (Riα RjβMαβ) (R-1ac (i)c) (R-1bd (j)d) = (Riα RjβMαβ) (R-1acδi,c) (R-1bdδj,d) = (Riα RjβMαβ) (R-1ai) (R-1bj) = (R-1ai Riα) (R-1bj Rjβ)Mαβ) = Mab Now how does the word "view" fit in here? The primed frame definitely has axis '(j). And in the primed frame, we "see" components M'ij which we can imagine to be diagonal with the right rotation. Principal Axis Transformation (ie, diagonalization of M) and a theorem. Imagine a rank-2 tensor M with components Mij in reference frame S having unit vectors (i). Imagine that R is an orthogonal transformation which brings M to diagonal form such that M' is diagonal. M' = RMRT // previous line in matrix notation M'ab = Raα RbβMαβ // contravariant rank 2 tensor In reference frame S' the tensor M has components M'ij and is diagonal, and this frame has basis vectors which are '(i) = R-1 (i). In reference frame S the tensor M has components Mij and is not diagonal, and this frame has basis vectors which are (i). What can be said about certain Mij elements? // R-1ij = '(i)j Mab = R-1acR-1bd M'cd = R-1acR-1bd (δc,dλc) = λc R-1acR-1bc = λc R-1acR-1bc = λc '(a)c '(b)c That is to say Mab = Σc λc '(a)c '(b)c Maa = Σc λc ( '(a)c)2 Imagine that the rotation R was selected such that the eigenvalues are in descending order λ1 ≥ λ2 .....≥ λn and we make no claim yet about signs. If λ1 > 0, then we may conclude that Maa = λ1 ( '(a)1)2 ≤ λ1 for all a = 1,2...n I once thought that in the case all λi are positive, you could show that all the largest Maa would be less than λ1 and the second largest would be smaller than λ2 and so on, but this is wrong as this counterexample shows: Here we have a real, square symmetric matrix Mij, and we find the three eigenvalues. This matrix is seen to be positive definite. And we find that all diagonal Mii ≤ λ1 (= 5) as predicted. But we do not find that the second largest Mii which here is 2 is less than λ2 which is .64. How unique is the R which does our principle axis transformation? I think the only freedom left is that you can use R or -R, and that must means reflecting all three of the primed unit vectors in S' space. Theorem: Suppose the symmetric matrix M is positive definite. We know that Maa = Σc λc ( '(a)c)2 and therefore all the diagonal elements Maa are positive. All we can say about these positive elements is that all of them are less than the largest eigenvalue λ1. Comments: Suppose you use this transformation on a galaxy of points as Ng does where Mij is the covariance matrix. Then when you go to the new frame S', you can arrange for M'ij to have descending positive eigenvalues and you know that M'11 will be larger than ALL the Mij and that is all you really know. In frame S', you have found axes such that the variances step down monotonically, and you have no correlation between pairs of directions. Ng Presentation To make the Ng notation work, we have to make two notational changes. (1) We need to make the replacement R → R-1 (same as R → RT) everywhere above. Among the results we find are V' = R-1 V // V is a rank 1 tensor (vector) with respect to R-1 M' = R-1MR = RTMR // M is a rank 2 tensor with respect to R-1 '(i) = R (i) // each (i) is a vector with respect to R (j) '(i) = '(i)j = Rji Now after making this change we have the '(i) being the columns of R R = ['(1), '(2), .... '(n)] . The form M' = R-1MR is the form one most often sees for a diagonalizing transformation, where the inverse is the matrix to the left of M. (2) We change the names of the basis vectors as follows (j) → j where j is a unit vector along the jth feature axis '(2) → (j) remove the prime. Notice that the components of vector j are (j)i and that (j)i = δj,i. Remember this just says that we have 1 = (1,0,0...0), 2 = (0,1,0...0), etc. Then the above equation set becomes V' = R-1 V // V is a rank 1 tensor (vector) with respect to R-1 M' = R-1MR = RTMR // M is a rank 2 tensor with respect to R-1 (i) = R i // each (i) is a vector with respect to R j (i) = (i)j = Rji R = [(1), (2), .... (n)] S = frame of reference having i as the unit vectors Su = frame of reference having (i) as the unit vectors For Ng, Mij is a positive definite covariance matrix. We rotate the basis vectors according to (i) = R i and in that new frame of reference which we might call Su we find that M has the diagonal form M'. The SVD. Ng calls upon the "singular value decomposition" which in the most general case and with a standard notation says M = USVH M = U S VH n x p n x n n x p p x p where we show the row x column sizes of the matrices. In general the elements of the matrices are complex numbers, and we have used H to indicate the Hermitian conjugate , as in VH = (VT)* or VHij = Vij* where * means complex conjugate. There are several notations for complex conjugate and for Hermitian conjugate and they conflict somewhat with one another complex conjugate use *; use overbar Hermitian conjugate use † ; use H ; use * ; use + In this decomposition, the square matrices U and V are both unitary, which means UH = U and VH = V. The non-square matrix S is all zeros except for some real numbers on its "diagonal" which starts in the upper left corner. These numbers are real and are called "the singular values" of matrix M and they are conventionally listed in monotonic decreasing order. If we write the U and V matrices in terms of their columns in this manner U = [(1), (2), .... (n)] V = [(1), (2), .... (n)] then the (i) are eigenvectors of MMH and the (i) are eigenvectors of MHM. Notice that both these matrices MMH and MHM are Hermitian (A=AH) and therefore have eigenvectors. For real numbers only, the condition is that A = AT for a matrix to have eigenvectors. The claim is that the non-zero singular values (on the diagonal of S) are the square roots of the non-zero eigenvalues of MMH and MHM. Our application of the SVD decomposition is much simpler than this general case. First of all, we have only real numbers so AH = AT for any of our matrices. Secondly, dimension p = n, so all matrices are square n x n. In this case, then, the SVD says M = U S VT where U and V are symmetric and S is diagonal. In Ng's application of SVD, we identify U = V = R where R is real orthogonal so R-1 = RT, and we identify S with M', our diagonal matrix, so the decomposition is then M = U S VT = R M' RT since M' = RTM R (see above) Since M is symmetric, it has some eigenvectors (i) and eigenvalues λi. Then we can verify the claim made above that the diagonal elements of S should be the square roots of the eigenvalues of MMT and MTM. Of course these two matrices are the same and are equal to M2 so: MMT (i) = MTM (i) = M2 (i) = M λi(i) = λi M(i) = λi λi(i) = λi2 (i) . If M is a positive definite matrix, then the descending elements of S = M' are the eigenvalues of M, λ1 ≥ λ2...... ≥ λn > 0 and the eigenvectors which are the columns of U = R are the corresponding eigenvectors. Main Point: We are trying to find the rotation matrix R which diagonalizes M into M' M' = RTM R such that the diagonal elements of M' are in descending order. By calling the SVD routine with M as argument and asking for three returned matrices [U, S, V] = svd (M) meaning M = U S VT = R M' RT when U is our desired matrix R, and S is our desired matrix M'. For the matrix M = Σ, the covariance matrix, we have [U, S, V] = svd (Σ) meaning Σ = U Σ' VT = R Σ' RT where Σ' is the covariance matrix observed in frame Su which is diagonal. Comment: Before diagonalization, a positive definite matrix M has all Mii > 0. After diagonalization the new diagonal M' matrix has diagonal elements M'ii = λi. These are all positive, and we know that the first one has the property that M'11 ≥ Mii for all i. Secondly, we know that M'ij = 0 when i and j are different, and this says that when viewed from frame Su, there is no correlation between any pair of different directions. Comment: Since the variances M'ii decrease as i increases, we want to consider dropping dimensions at the end of our list, not the beginning of the list. Reduction Idea: Ng says to do this for all the training samples x(i) z(i) = (Ureduced)T x(i) k k x n n [ Note conflict now, x(i) is a training sample, while (i) is an axis in feature space. ] Well here is what the means. We could start by doing this x'(i) = UT x(i) = R-1 x(i) and x'(i) is then our point x(i) transformed to x'-space, so to speak, according to V' = R-1 V in the Ng convention described above for R. The first k components of each x'(i) form a k-dimensional vector which Ng calls z(i) . Then we get exactly the above equation. This is the compressed feature data! Only the dimensions with the least variance were thrown out in this compression. We know that x(i) = U x'(i) If we replace U by Ureduced, we get some approximation to the original data x(i)approx = Ureduced x'(i) = Ureduced z(i) n n x n n n x k k How much "variance" was lost in our process. Here is the measure fraction lost = (1/m) Σi | x(i) - x(i)approx|2 / (1/m)Σi | x(i)|2 = (1/m) Σi | x(i) - x(i)approx|2 / Σj var(Xj) where the denominator is the sum of the variances over all the dimensions of feature space, and the numerator is the sum of the variance error in all dimensions caused by the imperfect reconstruction. Now how do we show that (1/m) Σi | x(i) - x(i)approx|2 / (1/m)Σi | x(i)|2 = 1 - ( Σi=1kλi/ Σi=1nλi) I haven't done much with variances with vectors inside, but see you can just add components. I do know that the trace of M is same as that of M', something I never used! M11 + M22 + .... + Mnn = λ1 + λ2 + ... + λn Well, here is a shot at it. x'(i) = UT x(i) = R-1 x(i) x(i)approx = Rr x'(i) = Rr R-1 x(i) where R = [(1), (2), .... (n)] Rr = [(1), (2), .... (k), 0,0,0..0] I think that RrRT = ((1)) ((1))T + ((2)) ((2))T + ... + ((k)) ((k))T Then we have x(i) - x(i)approx = x(i) - Rr R-1 x(i) = (1 - Rr R-1) x(i) = (1 - {((1)) ((1))T + ((2)) ((2))T + ... + ((k)) ((k))T }) x(i) = Σs=k+1n ((s)) ((s))T x(i) Then (1/m) Σi | x(i) - x(i)approx|2 = (1/m) Σi | (1 - Rr R-1) x(i)|2 = (1/m) Σi (1 - Rr R-1) x(i) (1 - Rr R-1) x(i) Back up. Go back to M = Mij (i) (j)T analogous to V = Vi (i) M = M'ij '(i) '(j)T analogous to V = V'i '(i) and write these in Ng notation as M = Σi,j=1nMij (i) (j)T M = Σi,j=1nM'ij (i) (j)T where (i) = R (i) Now I think he is setting the higher order u to zero, so we can say M ≈ Σi,j=1kM'ij (i) (j)T = Σi=1k M'ii(i) (i)T So with full M, we have total variances being Σs=1n λs and for approximate M we have Σs=1k λs, so this gives the notion that variance retained = Σs=1k λs / Σs=1n λs variance lost = 1 - Σs=1k λs / Σs=1n λs Now my only issue is how do you relate this to the other form. Well, I see now that Σs=1n λs = Σs=1n Mss since trace is preserved so you can regard either sum as the "total variance" as measured in either reference frame. So the question remaining is this: why does (1/m) Σi | x(i) - x(i)approx|2 represent the variance lost? It is reasonable, but I want a proof that it agrees with the λ stuff. Here is what I want to prove: (1/m) Σi | x(i) - x(i)approx|2 / Σs=1n λs = 1 - Σs=1k λs / Σs=1n λs or (1/m) Σi | x(i) - x(i)approx|2 = Σs=1n λs - Σs=1k λs = Σs=k+1n λs I did show above that x(i) - x(i)approx = Σs=k+1n ((s)) ((s))T x(i) but I don't see any nice way to use this. How about saying the variance is independent of rotated basis, so we have (1/m) Σi=1n | x(i) - x(i)approx|2 = (1/m) Σi=1n | x'(i) - x'(i)approx|2 = (1/m) Σi=1n Σs=1n(x's(i) - x'(i)approx,s)2 = (1/m) Σi=1n Σs=k+1n(x's(i) - x'(i)approx,s)2 = (1/m) Σi=1n Σs=k+1n(x's(i))2 = Σs=k+1n (1/m) Σi=1n(x's(i))2 I need to show that Σs=k+1n (1/m) Σi=1n(x's(i))2 = Σs=k+1n λs which means I need to show that (1/m) Σi=1n(x's(i))2 = λs But this is precisely the variance in the sth component direction. QED