# Issue with LinearAlgebra's `\`

I’m having some trouble with the `\` operator for matrices. Supposedly, for two matrices `A` and `B` it should give `x` satisfying `Ax = B`. However,

``````julia> A = [0 1 4; 1 3 5; 3 7 7]
3×3 Array{Int64,2}:
0  1  4
1  3  5
3  7  7

julia> B = [-5, -2, 6]
3-element Array{Int64,1}:
-5
-2
6

julia> x = A\B
3-element Array{Float64,1}:
5.254199565265581e15
-3.002399751580329e15
7.50599937895081e14

julia> A*x == B
false
``````

What is going on here?

1 Like

Have you checked what A*x actually gives?

notice your solution is at ~`e15` order… there may not be solution within Float64 range.

``````julia> A = [0 1 4; 1 3 5; 3 7 BigInt(7)]
3×3 Array{BigInt,2}:
0  1  4
1  3  5
3  7  7

julia> x = A\B
3-element Array{BigFloat,1}:
-6.7545385388434447330416407921734612914407491054956995689683590671282658956627e+76
3.8597363079105398474523661669562635951089994888546854679819194669304376546647e+76
-9.649340769776349618630915417390658987772498722136713669954798667326094136663e+75

julia> A*x
3-element Array{BigFloat,1}:
-5.0
-2.0
7.0
``````

maybe there’s a bug…

There are a couple of things going on here. First, the \ operator uses floating point arithmetic. So you can only expect A*x to be approximately equal to B. Secondly, your matrix appears to be singular, or numerically very close to being singular. In julia I get

``````julia> rank(A)
2
``````

In this case, the \ manages to produce a solution but it’s meaningless. Here’s an example with a full rank matrix, where you can expect a reasonably accurate solution:

``````julia> A = randn(3,3)
3×3 Array{Float64,2}:
0.296735   -0.708839   1.43307
-0.210209   -0.749724  -0.876588
-0.0215593   0.264439  -0.552682

julia> x = A\B
3-element Array{Float64,1}:
44.14624174273415
3.2040390103063894
-11.04520981290519

julia> A*x ≈ B
true

julia> A*x
3-element Array{Float64,1}:
-5.000000000000001
-1.9999999999999984
6.0

julia> A*x == B
false
``````

Note that in this case the solution is accurate but `A*x == B` still returns false because the solution is only accurate up to machine precision. It’s not the exact solution. I suggest reading up on numerical linear algebra to get a better sense of what’s going on here.

5 Likes

I knew that the system was inconsistent so I was expecting an error to be thrown instead of a nonsense value to be returned. For some other examples of inconsistent systems the operator throws a `SingularException(n)`, such as

``````julia> A = [1 1; 4 4]
2×2 Array{Int64,2}:
1  1
4  4

julia> B = [3, 10]
2-element Array{Int64,1}:
3
10

julia> A\B
ERROR: SingularException(2)
``````

However it was not throwing it for this one.

I was curious how other languages would handle the linear system.
Matlab gave a warning that the matrix is singular to working precision and returned `[NaN, Inf, -Inf]`.
Numpy with linalg.lstsq gave `[ 1.43722944, 2.05627706, -1.83549784]` but it also returns the rank with the fit (2 in this case). 1 Like

For the original matrix,

``````julia> lu(A)
LU{Float64,Array{Float64,2}}
L factor:
3×3 Array{Float64,2}:
1.0       0.0       0.0
0.0       1.0       0.0
0.333333  0.666667  1.0
U factor:
3×3 Array{Float64,2}:
3.0  7.0   7.0
0.0  1.0   4.0
0.0  0.0  -8.88178e-16
``````

which does not trigger the singular exception from LAPACK.

Generally you can’t rely on LAPACK catching ill-conditioning, it’s something you explicitly check for or design your algorithm to avoid.

2 Likes

The analogous thing in Julia is to use pivoted QR to do the minimum-norm least-square solution:

``````julia> qr(A, Val(true)) \ B
3-element Array{Float64,1}:
1.4372294372294387
2.056277056277058
-1.8354978354978373
``````

and from the `qr` factorization you can also see that it is rank deficient:

``````julia> qr(A, Val(true)).R
3×3 Array{Float64,2}:
-9.48683  -7.16783  -2.74064
0.0       2.76084   1.57762
0.0       0.0      -2.22045e-16
``````

(note the tiny third diagonal entry).

If you are doing linear algebra with matrices that might be ill-conditioned, you really need to think carefully about why that is happening in order to choose an appropriate course of action.

3 Likes