What's the best (in terms of computation time and numerical robustness) way to find the least squares solution of $$Ax = b$$ if $A$ is symmetric and positive semi-definite?
If $A$ were symmetric and positive definite, the Cholesky decomposition would seem to be the way to go to build the inverse. But that doesn't work for positive semidefinite cases.
Performing a spectral decomposition to get $A = Q \Lambda Q^T$, and then forming the pseudoinverse $A^+ = Q \Lambda^+ Q^T$ to get $x = A^+ b$ works, but I'm not convinced it's the most efficient.
According to MATLAB's
mldivide():If the Cholesky Decomposition doesn't work you should use LDL Decomposition.