1

I use ojAlgo to solve a system of linear equations. In one case I get a RecoverableCondition exception. Probably because matrix is ill-conditioned, the condition number is about 1e15.

I use ojAlgo to solve it as seen in the code below. It usually works, but not in this case.

Is there any other solver I could use for a symmetric indefinite (ill-conditioned) matrix?

The present failing size is 18x18 but later 1000x1000 might be needed. Since its part of a iterative algorithm the accuracy is not super important.

         SolverTask<Double> equationSolver = SolverTask.PRIMITIVE.make(KKT, rhs.negate());
         MatrixStore<Double> deltaX = null;
         try {
            deltaX = equationSolver.solve(KKT, rhs.negate());
         } catch (RecoverableCondition ex) {
            int i = 0;
         }

I tried to reproduce this in a self contained example but failed, because there it works. Maybe I do not get exactly the same matrix down to the last bit.

Magnus
  • 35
  • 5

1 Answers1

1

In your case, that method would use a Cholesky decomposition as the solver.

If here's a problem then try to pick another decomposition by instantiating a suitable alternative directly. An SVD can usually handle anything, but that would be very expensive. Perhaps QR can be ok.

QR<Double> qr = QR.PRIMITIVE.make(templateBody);
qr.decompose(body);
qr.getSolution(rhs,x);

This way you can reuse the decomposition instance as well as the solution vector x.

Another alternative is to precondition the body/KKT matrix. Perhaps add a small diagonal - just enough to make the Cholesky decomposition solvable.

if (!cholesky.isSolvable()) {
    // Fix that
}

Or perhaps try something in the org.ojalgo.matrix.task.iterative package.

apete
  • 1,250
  • 1
  • 10
  • 16
  • Thank you for your quick answer! I will try another decomposition, and/or preconditioning. According to the Templates book beginning of section 2.3.2 http://www.netlib.org/templates/templates.pdf CG might get into problem for symmetric indefinite systems. For Jacobi and Gauss–Seidel method convergence is only guaranteed for a limited class of matrices. – Magnus Oct 15 '21 at 09:18
  • I do not find MINRES (symmetric indefinite) or GMRES (general indefinite) in ojAlgo. Should they be added to ojAlgo in a identical way to CG? Maybe like http://netlib3.cs.utk.edu/templates/matlab/gmres.m https://web.stanford.edu/group/SOL/software/minres/ On the last page of the Templates book there is a flow chart for selecting an iterative method. It indicate CG might work "anyway". Cheers, Magnus – Magnus Oct 15 '21 at 09:18
  • That would be great contributions. – apete Oct 20 '21 at 05:51