JuliaLang / julia

The Julia Programming Language
https://julialang.org/
MIT License
45.68k stars 5.48k forks source link

Bugs in linear algebra. A * inv(A) != I, wrong determinants, wrong adjoints as well. #32845

Closed 0xDaksh closed 5 years ago

0xDaksh commented 5 years ago

The latest version of julia, has problems with linear algebra.

Versioninfo():

Julia Version 1.1.1
Commit 55e36cc308 (2019-05-16 04:10 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i5-8300H CPU @ 2.30GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.1 (ORCJIT, skylake)

Replication Script:

using LinearAlgebra

 A = [1 2 3; 4 5 6; 7 8 9]
 A * inv(A)

Out:

output

Wrong Determinants as well:

julia> det(A)
6.661338147750939e-16

Clearly the det(A) = 0

Wrong Adjoints as well:

julia> adjoint(A)
3×3 Adjoint{Int64,Array{Int64,2}}:
 1  4  7
 2  5  8
 3  6  9

adjoint of a

andreasnoack commented 5 years ago

As you point out, the matrix A is singular so I'm not sure why you'd expect A*inv(A) to be the identity. Ideally, it would give a SingularException and indeed that is what I get

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

julia> inv(A)
ERROR: SingularException(3)
Stacktrace:
 [1] checknonsingular at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/factorization.jl:12 [inlined]
 [2] #lu!#103(::Bool, ::Function, ::Array{Float64,2}, ::Val{true}) at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/lu.jl:41
 [3] #lu! at ./none:0 [inlined]
 [4] #lu#107 at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/lu.jl:142 [inlined]
 [5] lu at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/lu.jl:142 [inlined] (repeats 2 times)
 [6] inv(::Array{Int64,2}) at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/dense.jl:732
 [7] top-level scope at none:0

However, we round integer matrices to floating point when doing rational operations such as computing the inverse. This also explains why the determinant isn't zero. If you want to avoid floating points you can conveert to rationals

julia> Rational.(A)
3×3 Array{Rational{Int64},2}:
 1//1  2//1  3//1
 4//1  5//1  6//1
 7//1  8//1  9//1

If you inv that matrix then you should definitely get a SingularException.

Finally, it looks like you are expecting adjoint to return https://en.wikipedia.org/wiki/Adjugate_matrix but that is not what the documentation says

help?> adjoint
search: adjoint adjoint! Adjoint

  adjoint(A)

  Lazy adjoint (conjugate transposition) (also postfix '). Note that adjoint is applied recursively to elements.

  This operation is intended for linear algebra usage - for general data manipulation see permutedims.

Please see https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/index.html for more details about linear algebra in Julia or feel free to ask question on discourse.julialang.org.

0xDaksh commented 5 years ago

@andreasnoack yeah but this is only happening in the current stable release not in lts in which it throws the error, so I believe the expected behavior should be the singular exception which is not what is happening.

andreasnoack commented 5 years ago

so I believe the expected behavior should be the singular exception which is not what is happening.

You can't rely on the pivots being exactly zero when using floating point values. The tiniest rounding error could make the U factor non-singular. You have to use rational elements if you want to reliably detect exact singularity.