New answers tagged matrix
3
tch already gave you a good answer, but I have a suggestion for a simpler problem that you could address pretty easily and which should help you sanity check your results.
Let
$$M = \sum_i\lambda_iv_iv_i^\top$$
be the eigenvalue decomposition of $M$ and furthermore take the eigenvalues to be sorted in decreasing order: $\lambda_1 \ge \ldots \ge \lambda_n$.
...
6
This is a slightly modified version of my response on math.stackexchange. One standard approach to computing matrix functions times a vector $f(M)x$ or quadratic forms $x^Tf(M)x$ when $M$ is symmetric is via the Lanczos algorithm. Lanczos computes an orthonormal basis $Q_k = [q_1, \ldots, q_k]$ for Krylov subspace $\operatorname{span}(x,Ax,\ldots, A^{k-1}x)$ ...
3
For symmetric, more generally, normal, matrices there is no danger in diagonalizing a matrix to evaluate a matrix function, because all eigenvectors are perfectly conditioned. You should get answers with errors close to the machine epsilon times the condition number of the function evaluation problem. This is as accurate as you can get in floating point. ...
1
As mentioned, netlib BLAS is not at all optimized, but it is definetly the "refblas". Using IKML, ACML, OpenBLAS or "your vendor" BLAS, you are (somehow) assured, that the results of the operation of the optimized BLAS is equal to the "refblas" up to a known error. Take into care that: vendors (intel, amd, nvidia, ...) try hard ...
1
If you're ok with positive semi-definite matrices, then you can refer to this paper, which I think is what Frederico's answer also explains:
Higham NJ. Computing a nearest symmetric positive semidefinite matrix. Linear Algebra and its Applications. 1988 May;103(C):103-118.
A python implementation can be found in this answer, which also provides a link to a ...
Top 50 recent answers are included
