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?
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.
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).
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.
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.