Re: Matrix Decomposition

From: Gottfried Helms (helms_at_uni-kassel.de)
Date: 03/16/05 

Date: 15 Mar 2005 21:16:21 -0500

Am 15.03.05 22:15 schrieb Zdislav V. Kovarik:
> On Mon, 14 Mar 2005, Ben wrote:
>
>
>>Hi,
>>
>> We all know a matrix A can be decomposed into Q * R
>>
>>where Q is an orthogonal matrix and R is a triangular matrix.
>>
>>Any one know if we can find a decomposition Q * M for a matrix A,so
>>
>>A = Q * M
>>
>>where q is an orthogonal matrix and M is a positive definite symmetric matrix.
>>
>>Under what condition, if any, we can find such a decomposition??
>>
>>Your reply/comment will be hightly appreciated
>>
>>Ben
>
>
> Look up "polar decomposition". If A is a real square matrix, and "positive
> definite" admits the limit situation (some eigenvalues may be zero), there
> is no further condition. Uniqueness is not promised; it is guaranteed if A
> is invertible.
>
> One "overkill" method is: find the SVD: A = U * D * V^t, then
> Q = U * V^t, and M = V * D * V^t (t denotes transpose).
>
> One interesting property: If W runs through all orthogonal matrices
> then the trace of W^t * A is maximal for W=Q. Then, M = Q^t * A.
>
> Cheers, ZVK(Slavek)

Good information!

One other shot:
 poss 1) PC-decomposition can be done by pairwise rotation of columns,
         (Jacobi-Rotation) with a certain maximizing criterion. One can change
         that criterion to approximate a symmetric matrix instead
.
 poss 2) You do a procrustres rotation with the unit matrix as a target.

 poss 3) You build B=A*A and get the square-root with the matrix-analogon
          of the newton iteration.

Concerning the non-uniqueness of the solution: I prefer two extreme
solutions. The first where all entries in the diagonal are equal or
with least variance.
The second (and alternative) with maximum variance, ordered for high
absolute values descending from top left.

Gottfried Helms

;; Note I'm using T instead of Q usually (and here, sorry).
;===========================================
;MatMate-Listing vom:16.03.05 01:27:55
;============================================
      a :
             | 2.87 0.87 8.86 |
             | 0.48 6.21 7.90 |
             | 8.64 2.00 5.63 |

[1] t = gettrans(a',"ziel",einh(3))' //get the procrustes-rotation to Identity-Matrix
      t :
             | 0.39 -0.44 0.80 |
             | 0.91 0.09 -0.40 |
             | 0.11 0.89 0.44 |

[2] m = t*a // apply the rotation to A
      m :
             | 7.87 -0.80 4.53 |
             | -0.80 0.52 6.51 |
             | 4.53 6.51 10.47 |

// Newton-iteration to matrix-rooot
[23] b = a*a
[24] m = einh(3)

[25] m = (b*inv(m)+m)/2
[26] m = (b*inv(m)+m)/2
[27] m = (b*inv(m)+m)/2
[28] m = (b*inv(m)+m)/2
[29] m = (b*inv(m)+m)/2
[30] m = (b*inv(m)+m)/2
[31] m = (b*inv(m)+m)/2
[32] m = (b*inv(m)+m)/2
[33] m = (b*inv(m)+m)/2
[34] m = (b*inv(m)+m)/2
[35] m = (b*inv(m)+m)/2
[36] m = (b*inv(m)+m)/2
[37] m = (b*inv(m)+m)/2

      m :
             | 7.41 4.13 3.93 |
             | 4.13 8.84 2.45 |
             | 3.93 2.45 9.43 |

[38] t = m*inv(a) // What was the transformation ??

      t :
             | -0.02 -0.21 0.98 |
             | -0.46 0.87 0.18 |
             | 0.89 0.45 0.11 |

[39] chk = t*t' // T should be orthonormal

      chk :
             | 1.00 -0.00 -0.00 |
             | -0.00 1.00 0.00 |
             | -0.00 0.00 1.00 |

-------------------------------------