Hello,

I am trying to solve a linear system Ax = b using PartialPivLU and FullPivLU. For some reason, only PartialPivLU is capable of computing the solution.

I have already tried to adjust the threshold with a very small value, but no success. The maxPivot of the matrix is 3.791125e+18 and it seems A also has pivots in the order of 1.0e-2.

Below are the coefficient matrix A and the vector b. Why FullPivLU fails to solve such linear system?

Code: Select all

A =
-2.553787e+00 0 0 0 0 0 0 0 0 -1.801542e-02 1.783604e+03 0 9.367796e+02 0 0 0 0 0 0 0 0 0
0 -3.099768e+09 0 0 0 0 0 0 0 0 1.783604e+03 0 -9.367796e+02 0 0 0 0 0 0 0 0 0
0 0 -5.294874e+15 0 0 0 0 1 0 0 1.783604e+03 0 -9.367796e+02 0 0 0 0 0 0 0 0 0
0 0 0 -2.860755e+08 0 0 0 1 5.034438e-01 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 -4.795520e+08 0 0 1 2.517219e-01 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 -7.331359e+10 0 1 0 0 1.783604e+03 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 -3.791125e+18 1 0 1.801542e-02 0 0 -9.367796e+02 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 2.517219e-01 0 0 -3.155253e+03 9.367796e+02 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 1.783604e+03 0 -9.367796e+02 -4.357850e+09 0 0 0 2.517540e-01 0 0 0 0
0 0 0 0 0 0 0 0 2.517219e-01 0 0 0 0 0 -4.933290e-01 0 0 2.517540e-01 0 0 0 0
0 0 0 0 0 0 0 0 0 0 1.783604e+03 0 0 0 0 -1.213907e+03 0 2.517540e-01 0 0 0 0
0 0 0 0 0 0 0 0 0 1.801542e-02 0 0 -9.367796e+02 0 0 0 -3.014179e+10 2.517540e-01 0 0 0 0
0 0 0 0 0 0 0 0 0 1.801542e-02 0 0 -9.367796e+02 0 0 0 0 0 -1.766116e+10 0 0 0
0 0 0 0 0 0 0 0 0 1.801542e-02 0 0 0 0 0 0 0 0 0 -6.195980e+01 0 0
0 0 0 0 0 0 0 0 0 -1.801542e-02 1.783604e+03 0 9.367796e+02 0 0 0 0 0 0 0 -8.802000e-01 0
0 0 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 2 1 0 0 0 1 0 0 1 0 0 1 0 0 0 0 0 0 0
-1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 1 0 1 1 -1 0
1 1 1 0 0 1 0 0 0 0 1 0 0 1 0 1 0 0 0 0 1 0
1 -1 -1 0 0 0 -1 0 0 0 0 1 1 -1 0 0 -1 0 -1 0 1 0
0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1


Code: Select all

b =
1.776357e-15
3.552714e-15
1.065814e-14
1.665335e-14
2.664535e-15
4.440892e-15
-2.465583e-12
1.398881e-14
3.552714e-15
1.199041e-14
1.143530e-14
7.105427e-15
0
0
3.552714e-15
1.297455e-07
0
0
1.297455e-07
-1.297455e-07
0
0


Hi, according to double precision accuracy, your matrix is singular (rank 12), so the fact that partial piv LU gives you a better solution here is probably just luck.

The same matrix can be solved with GSL with no problems. This matrix is in fact a Jacobian matrix used during a Newton-Raphson iteration.

If I am not wrong, GSL performs a LU decomposition with partial pivoting. Perhaps this is why it succeeds in computing the linear system.

Is there a remedy that make it possible to solve the linear system Ax = b with FullPivLU?

forget what I said, there might be a threshold issue in FullPivLU, I added an entry in our bug tracker: http://eigen.tuxfamily.org/bz/show_bug.cgi?id=379

allanmulin wrote:I am trying to solve a linear system Ax = b using PartialPivLU and FullPivLU. For some reason, only PartialPivLU is capable of computing the solution.

Exactly what do you mean by "is capable of computing a solution"? How do you check this?

allanmulin wrote:I have already tried to adjust the threshold with a very small value, but no success. The maxPivot of the matrix is 3.791125e+18 and it seems A also has pivots in the order of 1.0e-2.

I thought that the threshold is not used by solve(). However, Gael seems to think otherwise. Can you please clarify: does changing the threshold make any difference to the result of solve() ?

jitseniesen wrote:Exactly what do you mean by "is capable of computing a solution"? How do you check this?

The solution I obtain with PartialPivLU matches the solution obtained by GSL using a LU decomposition with partial pivoting. The use of FullPivLU results in a solution vector x that fails ||Ax - b|| < epsilon for small enough epsilon (1.0e-6).

jitseniesen wrote:I thought that the threshold is not used by solve(). However, Gael seems to think otherwise. Can you please clarify: does changing the threshold make any difference to the result of solve() ?

No. Changing the threshold in FullPivLU yields the same solution vector x, which again fails the above checking.
I also tried to use a threshold = 0 as shown in the thread below, but with no success.
viewtopic.php?f=74&t=88127&p=159349&hilit=FullPivLU+octave#p159349