The determinant of your failing A is nonzero:
julia> det(A)
1.5543122344752192e-15
Of course, given its magnitude, that's probably just floating-point error
noise - but it makes me guess that the matrix inversion fails (i.e. doesn't
fail) for the same reason; somewhere along the calculations, enough
floating point error is accumulated to make the matrix non-singular (albeit
still very ill-conditioned, as seen by the monstrously large numbers in the
inverse...). I'm not confident enough to definitely rule out a bug
somewhere, but I'd consider it pretty likely that floating point arithmetic
is the root cause in this instance =)
// Tomas
On Monday, December 29, 2014 3:17:06 AM UTC+1, Ben Zeckel wrote:
>
>
> I am attempting to learn Julia and am experimenting with the matrix
> function inv to calculate the inverse of some matrices. Some are singular
> and some are nonsingular.
> One case (below) gives me the wrong result but it is a simple example so I
> must be misusing Julia (I found a similar issue with large numbers in the
> basic factorial function in the documentation returning 0 which came down
> to me not realizing the overflow behavior immediately)
>
> These attempts works as expected
>
> julia> A = [1 2 3; 0 2 2; 1 2 3]; inv(A)
> ERROR: SingularException(3)
> in inv at linalg/lu.jl:149
> in inv at linalg/dense.jl:328
>
> julia> A = [1 1 1; 0 2 3; 5 5 1]; inv(A)
> 3x3 Array{Float64,2}:
> 1.625 -0.5 -0.125
> -1.875 0.5 0.375
> 1.25 0.0 -0.25
>
> julia> inv(A) * A
> 3x3 Array{Float64,2}:
> 1.0 0.0 0.0
> 0.0 1.0 0.0
> 0.0 0.0 1.0
>
>
> This does not
>
> jjulia> A = [2 1 -1; 1 -2 -3; -3 -1 2]; inv(A)
> 3x3 Array{Float64,2}:
> -4.5036e15 -6.43371e14 -3.21686e15
> 4.5036e15 6.43371e14 3.21686e15
> -4.5036e15 -6.43371e14 -3.21686e15
>
> julia> A * inv(A)
> 3x3 Array{Float64,2}:
> 0.0 -0.125 -0.5
> 0.0 0.75 0.0
> 0.0 0.0 0.0
>
> Double checking manually, in julia with rref, and with wolframalpha.com
> shows the expected results on this singular matrix
>
> julia> A = [2 1 -1; 1 -2 -3; -3 -1 2]; rref([A eye(3)])
> 3x6 Array{Float64,2}:
> 1.0 0.0 -1.0 0.0 0.142857 -0.285714
> 0.0 1.0 1.0 0.0 -0.428571 -0.142857
> * 0.0 0.0 0.0 * 1.0 0.142857 0.714286
>
> inv {{2,1,-1}, {1,-2,-3}, {-3,-1,2}}
>
> http://www.wolframalpha.com/input/?i=inv+%7B%7B2%2C1%2C-1%7D%2C+%7B1%2C-2%2C-3%7D%2C+%7B-3%2C-1%2C2%7D%7D
>
> Thanks for any help,
>
> Ben
>
>
