RFC : isposdef, bkfact, eigvals! for Symmetric Matrix

[amitjamadagni]

bkfact and eigvals! give the following error

_       _ _(_)_     |  A fresh approach to technical computing
(_)     | (_) (_)    |  Documentation: http://docs.julialang.org
_ _   _| |_  __ _   |  Type "?help" for help.
| | | | | | |/ _` |  |
| | |_| | | | (_| |  |  Version 0.5.0-dev+3084 (2016-03-10 02:15 UTC)
_/ |\__'_|_|_|\__'_|  |  Commit c73feea (1 day old master)
|__/                   |  x86_64-linux-gnu

julia> a = sprand(5,5,0.3)
5x5 sparse matrix with 7 Float64 entries:
[5, 1]  =  0.328742
[1, 2]  =  0.070309
[4, 3]  =  0.49377
[5, 3]  =  0.362077
[2, 4]  =  0.974666
[3, 5]  =  0.666859
[4, 5]  =  0.278323

julia> a = a + a.'
5x5 sparse matrix with 12 Float64 entries:
[2, 1]  =  0.070309
[5, 1]  =  0.328742
[1, 2]  =  0.070309
[4, 2]  =  0.974666
[4, 3]  =  0.49377
[5, 3]  =  1.02894
[2, 4]  =  0.974666
[3, 4]  =  0.49377
[5, 4]  =  0.278323
[1, 5]  =  0.328742
[3, 5]  =  1.02894
[4, 5]  =  0.278323

julia> a = Symmetric(a)
5x5 Symmetric{Float64,SparseMatrixCSC{Float64,Int64}}:
0.0       0.070309  0.0      0.0       0.328742
0.070309  0.0       0.0      0.974666  0.0
0.0       0.0       0.0      0.49377   1.02894
0.0       0.974666  0.49377  0.0       0.278323
0.328742  0.0       1.02894  0.278323  0.0

julia> bkfact(a)
ERROR: MethodError: no method matching bkfact(::SparseMatrixCSC{Float64,Int64}, ::Symbol, ::Bool)
Closest candidates are:
bkfact{T<:Union{Complex{Float32},Complex{Float64},Float32,Float64}}(::Union{DenseArray{T<:Union{Complex{Float32},Complex{Float64},Float32,Float64},2},SubArray{T<:Union{Complex{Float32},Complex{Float64},Float32,Float64},2,A<:DenseArray{T,N},I<:Tuple{Vararg{Union{Base.AbstractCartesianIndex{N},Base.NoSlice,Colon,Int64,Range{Int64}}}},L}}, ::Symbol, ::Bool)
bkfact{T<:Union{Complex{Float32},Complex{Float64},Float32,Float64}}(::Union{DenseArray{T<:Union{Complex{Float32},Complex{Float64},Float32,Float64},2},SubArray{T<:Union{Complex{Float32},Complex{Float64},Float32,Float64},2,A<:DenseArray{T,N},I<:Tuple{Vararg{Union{Base.AbstractCartesianIndex{N},Base.NoSlice,Colon,Int64,Range{Int64}}}},L}}, ::Symbol, ::Bool, ::Bool)
bkfact{T}(::Union{DenseArray{T,2},SubArray{T,2,A<:DenseArray{T,N},I<:Tuple{Vararg{Union{Base.AbstractCartesianIndex{N},Base.NoSlice,Colon,Int64,Range{Int64}}}},L}}, ::Symbol, ::Bool)
...
[inlined code] from ./expr.jl:8
in bkfact(::Symmetric{Float64,SparseMatrixCSC{Float64,Int64}}) at ./linalg/symmetric.jl:141
in eval(::Module, ::Any) at ./boot.jl:267

julia> det(a)
ERROR: MethodError: no method matching bkfact(::SparseMatrixCSC{Float64,Int64}, ::Symbol, ::Bool)
Closest candidates are:
bkfact{T<:Union{Complex{Float32},Complex{Float64},Float32,Float64}}(::Union{DenseArray{T<:Union{Complex{Float32},Complex{Float64},Float32,Float64},2},SubArray{T<:Union{Complex{Float32},Complex{Float64},Float32,Float64},2,A<:DenseArray{T,N},I<:Tuple{Vararg{Union{Base.AbstractCartesianIndex{N},Base.NoSlice,Colon,Int64,Range{Int64}}}},L}}, ::Symbol, ::Bool)
bkfact{T<:Union{Complex{Float32},Complex{Float64},Float32,Float64}}(::Union{DenseArray{T<:Union{Complex{Float32},Complex{Float64},Float32,Float64},2},SubArray{T<:Union{Complex{Float32},Complex{Float64},Float32,Float64},2,A<:DenseArray{T,N},I<:Tuple{Vararg{Union{Base.AbstractCartesianIndex{N},Base.NoSlice,Colon,Int64,Range{Int64}}}},L}}, ::Symbol, ::Bool, ::Bool)
bkfact{T}(::Union{DenseArray{T,2},SubArray{T,2,A<:DenseArray{T,N},I<:Tuple{Vararg{Union{Base.AbstractCartesianIndex{N},Base.NoSlice,Colon,Int64,Range{Int64}}}},L}}, ::Symbol, ::Bool)
...
[inlined code] from ./expr.jl:8
in det(::Symmetric{Float64,SparseMatrixCSC{Float64,Int64}}) at ./linalg/symmetric.jl:146
in eval(::Module, ::Any) at ./boot.jl:267

julia> eigvals!(a)
ERROR: MethodError: no method matching eigvals!(::Symmetric{Float64,SparseMatrixCSC{Float64,Int64}})
in eval(::Module, ::Any) at ./boot.jl:267

The PR attempts to fix this :

julia> a = sprand(5,5,0.3)
5x5 sparse matrix with 8 Float64 entries:
    [1, 1]  =  0.978011
    [4, 1]  =  0.942458
    [1, 2]  =  0.370794
    [4, 2]  =  0.843244
    [5, 2]  =  0.972102
    [5, 3]  =  0.0898259
    [1, 4]  =  0.876004
    [2, 4]  =  0.0214496

julia> a = a + a.'
5x5 sparse matrix with 11 Float64 entries:
    [1, 1]  =  1.95602
    [2, 1]  =  0.370794
    [4, 1]  =  1.81846
    [1, 2]  =  0.370794
    [4, 2]  =  0.864694
    [5, 2]  =  0.972102
    [5, 3]  =  0.0898259
    [1, 4]  =  1.81846
    [2, 4]  =  0.864694
    [2, 5]  =  0.972102
    [3, 5]  =  0.0898259

julia> a = Symmetric(a)
5x5 Symmetric{Float64,SparseMatrixCSC{Float64,Int64}}:
 1.95602   0.370794  0.0        1.81846   0.0
 0.370794  0.0       0.0        0.864694  0.972102
 0.0       0.0       0.0        0.0       0.0898259
 1.81846   0.864694  0.0        0.0       0.0
 0.0       0.972102  0.0898259  0.0       0.0

julia> bkfact(a)
Base.LinAlg.BunchKaufman{Float64,Array{Float64,2}}(5x5 Array{Float64,2}:
 0.396454  0.428815  -22.7589     -0.0         0.381435
 0.370794  0.0        -0.079901   -0.0         0.889509
 0.0       0.0         0.0         0.0924038  -0.0
 1.81846   0.864694    0.0         0.0         0.972102
 0.0       0.972102    0.0898259   0.0         0.0     ,[1,-2,-2,-2,-2],'U',true,false,0)

julia> det(a)
0.002391777723342779

julia> eigvals!(a)
5-element Array{Float64,1}:
 -1.49922
 -0.707761
  0.000756814
  0.918071
  3.24417

There was a reference to the above features of Symmetric matrices in the following PR's
JuliaLang/LinearAlgebra.jl#253, JuliaLang/LinearAlgebra.jl#298,
#13166, JuliaLang/LinearAlgebra.jl#252

I hope this is going in the right direction and also I hope to work on the tests given that this is going in the expected direction. Hoping to hear from the community.

[andreasnoack]

Hi @amitjamadagni

Thanks for the PR. There are several reasons why you can get a MethodError in Julia. Sometimes there is a bug, sometimes a method hasn't been implemented yet and sometimes it doesn't necessarily makes sense to have the specific method defined. I'll add specific comments to your specific changes.

[andreasnoack]

Since Symmetric is a symmetric view of the matrix it is not correct to call isposdef on the data field. In many cases data will be a triangular matrix. Instead, you should probably expand this signature to allow for Symmetric and Hermitian sparse matrices.

[amitjamadagni]

I tried allowing more types to the signature but this did not work

function isposdef{Tv<:VTypes}(A::Union{SparseMatrixCSC{Tv,SuiteSparse_long}, Symmetric})
    if !ishermitian(A)
        return false
    end
    try
        f = cholfact(A)
    catch e
        isa(e, LinAlg.PosDefException) || rethrow(e)
        return false
    end
    true
end

and also there is a extra check ishermitian which is not required. I have pushed a commit with changes, I am not sure if this is the optimal way of doing things to add such a feature but any comment would be helpful. The reference for this is the following : understanding of this wiki article, which gives a iff condition for a matrix being positive definite in relation to Cholesky decomposition.

[andreasnoack]

Densifying the sparse matrices is not the right approach so I'll close.