On Mon, Feb 2, 2009 at 6:03 AM, Jitse Niesen
<jitse@xxxxxxxxxxxxxxxxx> wrote:
On Mon, 2 Feb 2009, Keir Mierle wrote:
I implemented a simple pivoting strategy for LDLT, but am not convinced it's better.
It looks for the largest diagonal element during each iteration, then pivots it to become the next processed row. Attached is a simple test that compares the old LDLT vs the pivoting one; it shows that the results are not fantastic. I find that there is consistently 30% or more residual in the pivoting version.
When I try it out on 100 random matrices of size 100 x 100, I get:
Average residual: nopiv 4.45718e-14 piv 4.07206e-14
Std.dev.: nopiv 1.43516e-15 piv 1.01811e-15
Minimum: nopiv 4.09651e-14 piv 3.77983e-14
Maximum: nopiv 4.83319e-14 piv 4.28309e-14
So pivoting is slightly better, but the residual matrix LDL^T - A is even without pivoting so small that it's not clear whether it's worth the effort for random matrices. However, random matrices often work surprisingly well.
I attach my code, which is the same as Keir but with some parts commented out that I did not understand or found unnecessary (please let me know if I commented out any essential parts).
Here is one example in which pivoting does vastly better:
A << 1e-5*MatrixXd::Identity(s/2,s/2), MatrixXd::Identity(s/2,s/2),
MatrixXd::Identity(s/2,s/2), 1e+5*MatrixXd::Identity(s/2,s/2);
A += 1e-5*MatrixXd::Random(s,s);
A *= A.transpose();
// and remove all the bits with SelfAdjointEigenSolver
I did some more experiments, and I'm convinced. Not only does the example you give here fail badly with nopivot, the nopivot version is not rank-revealing in the semidefinite case. What's left is the API decision; should I make pivoting a template option (enabled by default)? Since I doubt many LDLT objects get instantiated this is probably ok.
However, now the matrix L is not meaningful without the permutation matrix. This is not an issue for solving, but could be for users of .matrixL().. On the other hand, this is currently the case for LU.
Thoughts?
Keir
p.s. the blocked LDLT which can deal with indefinite matrices is a different algorithm (D is block-diagonal with 1x1 and 2x2 blocks), and I am not targetting it because I only need to solve semidefinite systems stably.
Incidentally, the documentation for setRandom() can be improved. It should mention what distribution is used and also something about setting the seed..
Cheers,
Jitse