Creating an equation library (Updated 03-FEB-2017)

[Han]

(01-15-2015 12:50 AM)Snorre Wrote:  Hello,
Not an answer to why complex elements in the first matrix: just observed that the second vector element is nearly zero, so the complex elements have no impact on L1(1)*diag(L1(2))*TRN(L1(3)).
Seems like a floating point rounding issue in the imaginary part calculation of Σ which translates to some really bad error in U, so that it isn't unitary anymore (L1(1)*trn(L1(1))=trn(L1(1))*L1(1)=identity doesn't hold).

Maybe the algorithm uses 1/s_i * A * v_i where s_i is the i-th eigenvalue and v_i is the corresponding eigenvector (perhaps this is why the corresponding column of U is all complex -- but I can't seem to make the results match this...)

Another way to get non-diagonal matrices is to compute the orthonormal eigenvectors of A^T * A and of A * A^T, for which the Prime does not seem to introduce complex values. The problem is that for non-square matrices additional steps may be needed to ensure that the eigenvectors are orthogonal (usually Gram-Shmidt takes care of this).

Still, that's a lot of additional work that should have been automatically handled by SVD... :-(