That probably confirms my hypothesis: when you move down to Float32, the
rounding errors from the arithmetic are truncated away, so that the
singularity of the matrix is preserved.
// T
On Monday, December 29, 2014 4:42:01 AM UTC+1, Ben Zeckel wrote:
>
> Interesting. if I move down to Float32 (not the default on my 64-bit
> version), I a better result.
>
> julia> versioninfo()
> Julia Version 0.3.4
> Commit 3392026* (2014-12-26 10:42 UTC)
> Platform Info:
> System: Windows (x86_64-w64-mingw32)
> CPU: Intel(R) Core(TM) i5-4200U CPU @ 1.60GHz
> WORD_SIZE: 64
> BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Haswell)
> LAPACK: libopenblas
> LIBM: libopenlibm
> LLVM: libLLVM-3.3
>
> julia> A = [2 1 -1; 1 -2 -3; -3 -1 2];
>
> julia> A
> 3x3 Array{Float64,2}:
> 2.0 1.0 -1.0
> 1.0 -2.0 -3.0
> -3.0 -1.0 2.0
>
> julia> det(A)
> 1.5543122344752192e-15
>
> julia> eps(Float64)
> 2.220446049250313e-16
>
> julia> nextfloat(0.0)
> 5.0e-324
>
> julia> A = float32(A)
> 3x3 Array{Float32,2}:
> 2.0 1.0 -1.0
> 1.0 -2.0 -3.0
> -3.0 -1.0 2.0
>
> julia> det(A)
> 0.0f0
>
> julia> inv(A)
> ERROR: SingularException(3)
> in inv at linalg/lu.jl:149
> in inv at linalg/dense.jl:328
>
> julia>
>
> On Sunday, December 28, 2014 10:05:15 PM UTC-5, Tomas Lycken wrote:
>>
>> 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
>>>
>>>
