# Issue with LinearAlgebra's \`\\\`

**URL:** <https://discourse.julialang.org/t/issue-with-linearalgebras/28901>\
**Category:** General Usage\
**Created:** [September 18, 2019, 5:01pm UTC](https://discourse.julialang.org/t/issue-with-linearalgebras/28901 "2019-09-18T17:01:31Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![srb](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/srb/32/6812_2.png) [@srb](https://discourse.julialang.org/u/srb)\
**Post date:** [September 18, 2019, 5:01pm UTC](https://discourse.julialang.org/t/issue-with-linearalgebras/28901/1 "2019-09-18T17:01:31Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [September 18, 2019, 5:05pm UTC](https://discourse.julialang.org/t/issue-with-linearalgebras/28901/2 "2019-09-18T17:05:32Z")

</div>

Have you checked what A\*x actually gives?

---

<div class="post-metadata">

**Author:** ![jling](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jling/32/212909_2.png) [@jling](https://discourse.julialang.org/u/jling)\
**Post date:** [September 18, 2019, 5:10pm UTC](https://discourse.julialang.org/t/issue-with-linearalgebras/28901/3 "2019-09-18T17:10:37Z")

</div>

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

```julia
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…

---

<div class="post-metadata">

**Author:** ![Pbellive](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pbellive/32/3604_2.png) [@Pbellive](https://discourse.julialang.org/u/Pbellive)\
**Post date:** [September 18, 2019, 5:14pm UTC](https://discourse.julialang.org/t/issue-with-linearalgebras/28901/4 "2019-09-18T17:14:15Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![srb](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/srb/32/6812_2.png) [@srb](https://discourse.julialang.org/u/srb)\
**Post date:** [September 18, 2019, 5:20pm UTC](https://discourse.julialang.org/t/issue-with-linearalgebras/28901/5 "2019-09-18T17:20:34Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![zaki](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zaki/32/10300_2.png) [@zaki](https://discourse.julialang.org/u/zaki)\
**Post date:** [September 18, 2019, 5:25pm UTC](https://discourse.julialang.org/t/issue-with-linearalgebras/28901/6 "2019-09-18T17:25:45Z")

</div>

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).  
🤷‍♂️

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [September 18, 2019, 5:34pm UTC](https://discourse.julialang.org/t/issue-with-linearalgebras/28901/7 "2019-09-18T17:34:06Z")

</div>

For the original matrix,

```julia
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.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [September 18, 2019, 5:46pm UTC](https://discourse.julialang.org/t/issue-with-linearalgebras/28901/8 "2019-09-18T17:46:59Z")

</div>

> [@zaki](#):
>
> Numpy with linalg.lstsq gave

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

```julia
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
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.