Check the solution based on the factorization Matrices of KLU

General Usage
1

I am trying to play with KLU matrices based on KLU.jl. However, I am having false instead of true when I am comparing Q*((P*A*Q)\P*b) == A\b. Could you please correct me given that the below A and b are taken from my real code?

A = sparse([0.0137346 0 0 0 0 0 0 0 0 -0.00123457 0 0 0 0 0;
            0 0.0137346 0 0 0 0 0 0 0 0 -0.00123457 0 0 0 0;
            0 0 0.0137346 0 0 0 0 0 0 0 0 -0.00123457 0 0 0;
            0 0  0 1.0 0 0 -1.0 0 0 0 0 0 1.0 0 0;
            0 0 0 0 1.0 0 0 -1.0 0 0 0 0 0 1.0 0;
            0 0 0 0 0 1.0 0 0 -1.0 0 0 0 0 0 1.0;
            0 0 0 -1.0 0 0 1.00302 -0.000705766 -0.000705766 0 0 0 0 0 0;
            0 0 0 0 -1.0 0 -0.000705766 1.00302 -0.000705766 0 0 0 0 0 0;
            0 0 0 0 0 -1.0 -0.000705766 -0.000705766 1.00302 0 0 0 0 0 0;
            -0.00123457 0 0 0 0 0 0 0 0 0.00425693 -0.000705766 -0.000705766 0 0 0;
            0 -0.00123457 0 0 0 0 0 0 0 -0.000705766 0.00425693 -0.000705766 0 0 0;
            0 0 -0.00123457 0 0 0 0 0 0 -0.000705766 -0.000705766 0.00425693 0 0 0;
            0 0 0 1.0 0 0 0 0 0 0 0 0 0 0 0;
            0 0 0 0 1.0 0 0 0 0 0 0 0 0 0 0;
            0 0 0 0 0 1.0 0 0 0 0 0 0 0 0 0]);

b = [-1.7632749965554524,4.5258858371825, -2.7626108406269907, 0.0, 0.0, 0.0, -84.69238104647695, 334.2439006136896, -249.55151956721267, 14.008262358676006, -19.104993786251043, 5.096731427574866,
     16329.20635908069, -8031.318664614193, -8297.88769446649];

factor = klu(A);

julia> (A[factor.p,factor.q]\b[factor.p])[factor.q] == A\b  # which is Q*((P*A*Q)\P*b) == A\b
false
2

This is just floating point error. Try with isapprox

1 Like
3

Thanks for your reply. It gives false as well. Actually, if we look at the resulting matrices, we can see that the order is different, in which it should be the same.

julia> isapprox((A[factor.p,factor.q] \ b[factor.p])[factor.q],A\b)
false

julia> A \b
15-element Vector{Float64}:
   128.21179193095762
   -16.895499498950347
  -111.31629243200852
 16329.20635908069
 -8031.318664614193
 -8297.88769446649
 16184.2153786398
 -7668.503713593524
 -8515.71166504627
  2854.607413115808
 -3853.9238476560936
   999.3164345402257
  -144.9909804408917
   362.8149510206699
  -217.82397057978096

julia> (A[factor.p,factor.q] \ b[factor.p])[factor.q]
15-element Vector{Float64}:
   128.21179193095762
   -16.895499498950347
  -111.31629243200852
 -7668.503713593524
 -8515.71166504627
 16184.2153786398
 16329.20635908069
 -8031.318664614193
 -8297.88769446649
  -144.9909804408917
   362.8149510206699
  -217.82397057978096
 -3853.9238476560936
   999.3164345402257
  2854.607413115808
5

Sorting the vectors before comparing might solve this issue.

6

Thanks for your feedback.
Based on which factor should I sort? Could you please guide me here cause I am not able to figure it out?

7

(A[factor.p,factor.q]\b[factor.p])[invperm(factor.q)] ≈ A\b
returns true (used the operator version of isapprox).
But I think you would want to use factor \ b, as this uses the specialized \ solve operator.

2 Likes
8

Got it.
Thank you very much for your valuable guide.