Re: [eigen] Pivoting for LDLT

```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.
```
That surprises me.

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

```
Incidentally, the documentation for setRandom() can be improved. It should mention what distribution is used and also something about setting the seed.
```
Cheers,
Jitse```
```#define EIGEN_DEFAULT_IO_FORMAT IOFormat(4, AlignCols, "  ", "\n", "", "", "", "")

#include <Eigen/Core>
#include <Eigen/Cholesky>
#include <Eigen/LU>
#include <Eigen/QR>
#include <Eigen/SVD>
#include <Eigen/Array>

using namespace Eigen;
using namespace std;

template<typename MatrixType>
void computeLDL_nopivot(const MatrixType& a, MatrixType &m_matrix)
{
assert(a.rows()==a.cols());
const int size = a.rows();
m_matrix.resize(size, size);
bool m_isPositiveDefinite = true;
typedef typename MatrixType::RealScalar RealScalar;
typedef typename MatrixType::Scalar Scalar;
const RealScalar eps = precision<Scalar>();

if (size<=1) {
m_matrix = a;
return;
}

Matrix<Scalar,MatrixType::RowsAtCompileTime,1> _temporary(size);

m_matrix.row(0) = a.row(0).conjugate();
m_matrix.col(0).end(size-1) = m_matrix.row(0).end(size-1) / m_matrix.coeff(0,0);
for (int j = 1; j < size; ++j)
{
RealScalar tmp = ei_real(a.coeff(j,j) - (m_matrix.row(j).start(j) * m_matrix.col(j).start(j).conjugate()).coeff(0,0));
m_matrix.coeffRef(j,j) = tmp;

if (0 && tmp < eps) {
m_isPositiveDefinite = false;
cout << "dying; with no pivot tmp = " << tmp << " rank = " << (j+1) << endl;
return;
}

int endSize = size-j-1;
if (endSize > 0) {
_temporary.end(endSize) = ( m_matrix.block(j+1,0, endSize, j)
* m_matrix.col(j).start(j).conjugate() ).lazy();

m_matrix.row(j).end(endSize) = a.row(j).end(endSize).conjugate()
- _temporary.end(endSize).transpose();

m_matrix.col(j).end(endSize) = m_matrix.row(j).end(endSize) / tmp;
}
}
}

// Computes a (PL)D(PL)^T decomposition. m_matrix holds the result P*L
// afterwords.
template<typename MatrixType>
void computeLDL_pivot(const MatrixType& a, MatrixType &m_matrix)
{
assert(a.rows()==a.cols());
const int size = a.rows();
m_matrix.resize(size, size);
bool m_isPositiveDefinite = true;
typedef typename MatrixType::RealScalar RealScalar;
typedef typename MatrixType::Scalar Scalar;
const RealScalar eps = precision<Scalar>();

if (size<=1) {
m_matrix = a;
return;
}

m_matrix = a; //k

typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntVectorType;
IntVectorType transpositions(size);
RealScalar biggest = RealScalar(0), biggest_in_corner;

Matrix<Scalar,MatrixType::RowsAtCompileTime,1> _temporary(size);

for (int j = 0; j < size; ++j)
{
// Find largest element diagonal and pivot it upward for processing next.
int row_of_biggest_in_corner, col_of_biggest_in_corner;
biggest_in_corner = m_matrix.diagonal().end(size-j).cwise().abs()
.maxCoeff(&row_of_biggest_in_corner,
&col_of_biggest_in_corner);
ei_assert(col_of_biggest_in_corner == 0); // XXX

if(j == 0) biggest = biggest_in_corner;

// Finish early if the matrix is not full rank.
if(0 && ei_isMuchSmallerThan(biggest_in_corner, biggest))
{
m_isPositiveDefinite = false;
for(int i = j; i < size; i++) transpositions.coeffRef(i) = i;
cout << "dying; not full rank. biggest in corner is much smaller than biggest";
break;
}

row_of_biggest_in_corner += j;
transpositions.coeffRef(j) = row_of_biggest_in_corner;
if(j != row_of_biggest_in_corner)
{
m_matrix.row(j).swap(m_matrix.row(row_of_biggest_in_corner));
m_matrix.col(j).swap(m_matrix.col(row_of_biggest_in_corner));
}

if (j == 0) {
m_matrix.row(0) = m_matrix.row(0).conjugate();
m_matrix.col(0).end(size-1) = m_matrix.row(0).end(size-1) / m_matrix.coeff(0,0);
continue;
}

RealScalar tmp = ei_real(m_matrix.coeff(j,j) - (m_matrix.row(j).start(j) * m_matrix.col(j).start(j).conjugate()).coeff(0,0));
m_matrix.coeffRef(j,j) = tmp;

if (0 && tmp < eps)
{
m_isPositiveDefinite = false;
for(int i = j; i < size; i++) transpositions.coeffRef(i) = i;
cout << "DYING; tmp = " << tmp << " rank = " << (j+1) << endl;
break;
}

int endSize = size-j-1;
if (endSize > 0) {
_temporary.end(endSize) = ( m_matrix.block(j+1,0, endSize, j)
* m_matrix.col(j).start(j).conjugate() ).lazy();

m_matrix.row(j).end(endSize) = m_matrix.row(j).end(endSize).conjugate()  //kk
- _temporary.end(endSize).transpose();

m_matrix.col(j).end(endSize) = m_matrix.row(j).end(endSize) / tmp;
}
}

// cout << "perms: " << transpositions.transpose() << endl;

IntVectorType m_p(size);
for(int k = 0; k < size; ++k) m_p.coeffRef(k) = k;
for(int k = size-1; k >= 0; --k) {
std::swap(m_p.coeffRef(k), m_p.coeffRef(transpositions.coeff(k)));
}
// cout << "m_p: " << m_p.transpose() << endl;

// For now, make a permuted L matrix so the caller can use it as A = LL^T,
// even though it's really PL, not L, that they are getting.
_temporary = m_matrix.diagonal();
MatrixType t_matrix = m_matrix.template part<Eigen::UnitLowerTriangular>();
m_matrix = t_matrix;
m_matrix.diagonal() = _temporary;
for (int j = size-1; j >= 0; --j) {
m_matrix.row(j).swap(m_matrix.row(transpositions.coeffRef(j)));
m_matrix.col(j).swap(m_matrix.col(transpositions.coeffRef(j)));
}

}

template<int s, typename T>
void tryDecomp(const T &A, double *res_nopiv, double *res_piv) {
// cout << "rank =\n" << A.lu().rank() << endl;
// cout << "Ap=\n" << A << endl;
// cout << "eigenvalues:\n" << (A.template part<SelfAdjointBit>().eigenvalues()) << endl;

Matrix<double, s, s> ldlt;
Matrix<double, s, s> L, D, res, X;

// No pivoting; original algorithm.
computeLDL_nopivot(A, ldlt);
L = ldlt.template part<Eigen::UnitLowerTriangular>();
D = ldlt.diagonal().asDiagonal();
res = L*D*L.transpose() - A;
// cout << "nopiv res = \n" << (res.cwise().abs().cwise() >= 1e-10) << endl;
// cout << "nopiv res = \n" << res << endl;
// cout << "nopiv res = " << sqrt(res.cwise().square().sum()) << endl;
*res_nopiv = sqrt(res.cwise().square().sum());

// With pivoting; doesn't work as well.
computeLDL_pivot(A, ldlt);
L = ldlt; // NOT lower triangular anymore (this is really PL).
L.diagonal().fill(1);
D = ldlt.diagonal().asDiagonal();
res = L*D*L.transpose() - A;
// cout << "piv res = \n" << (res.cwise().abs().cwise() >= 1e-10) << endl;
// cout << "piv res = \n" << res << endl;
// cout << "piv res = " << sqrt(res.cwise().square().sum()) << endl;
*res_piv = sqrt(res.cwise().square().sum());
}

int main() {
const int s = 100;
const int experiments = 100;

Matrix<double, s, s> A;
Matrix<double, s, 1> b;
Matrix<double, experiments, 1> res_nopiv, res_piv;

for (int i=0; i<experiments; i++) {
A.setRandom();
// A = A.cwise().abs();
A *= A.transpose();
// A /= 2;
// cout << "A=\n" << A << endl;
// cout << "eigenvalues orig:\n" << A.part<SelfAdjointBit>().eigenvalues() << endl;

// Force A to be badly conditioned, but p.s.d. No pivoting consistently works
// better.
b = eis.eigenvalues().cwise().abs();
b(0) *= 1e-9;
b(1) *= 1e-11;
b(4) *= 1e-11;
b(8) *= 1e-14;
A = eis.eigenvectors() * b.asDiagonal() * eis.eigenvectors().transpose();
// A(4, 4) = 0.0;
tryDecomp<s>(A, &res_nopiv[i], &res_piv[i]);

// N*N^t is on order 1e-60; residual error from decomposition is indeed
// 1e-75, which is about machine precision compared to 1e-60.
// Matrix<double, s, s> N;
// N = A.part<Eigen::UnitLowerTriangular>();
// N *= 1e-30;
// cout << "N \n" << N << endl;
// cout << "\n\n\n";
// cout << "This should work well with error about machine precision compared to the A matrix.";
// cout << "Note that the tmp < eps check in trunk LDLT prevents this from solving.";
// A = N * N.transpose();
// tryDecomp<s>(A);
}
cout << "Average residual:  nopiv " << res_nopiv.sum() / experiments;
cout << "    piv " << res_piv.sum() / experiments << endl;
cout << "Std.dev.:          nopiv " << sqrt(res_nopiv.cwise().square().sum() / experiments - pow(res_nopiv.sum() / experiments, 2));
cout << "    piv " << sqrt(res_piv.cwise().square().sum() / experiments - pow(res_piv.sum() / experiments, 2)) << endl;
cout << "Minimum:           nopiv " << res_nopiv.minCoeff();
cout << "    piv " << res_piv.minCoeff() << endl;
cout << "Maximum:           nopiv " << res_nopiv.maxCoeff();
cout << "    piv " << res_piv.maxCoeff() << endl;
}
```

 Mail converted by MHonArc 2.6.19+ http://listengine.tuxfamily.org/