From: Ian McCulloch (ianmcc_at_[hidden])
Date: 2007-05-03 17:22:53
Stjepan Rajko wrote:
> Philippe,
>
> The following code snippets might be helpful - I used singular value
> decomposition using lapack bindings to get eigenstuff. I don't know
> if it's the same thing as eig() in your case.
>
> ublas::matrix<double > A(rows, cols);
>
> ublas::matrix<double > ATA;
> ATA = prod(trans(A),A);
>
> ublas::matrix<double, boost::numeric::ublas::column_major > U(cols,
> cols), Vt(cols,cols);
> ublas::vector<double> S(cols);
> boost::numeric::bindings::lapack::gesvd(ATA, S, U, Vt);
>
> cout << U << endl << endl;
> cout << Vt << endl << endl;
> cout << S << endl << endl;
It would be easier and more accurate to calculate the singular value
decomposition of A itself, why bother with prod(trans(A),A) ? If you
really do want to diagonalize prod(trans(A),A), a symmetric (and
non-negative) matrix, the SVD is equivalent to an ordinary eigensolver and
you would be better off using syev or syevr.
Cheers,
Ian