Confused `spdiagm` with matlab function

Specific Domains Numerics
1

I’m having issues running cholfact on a sparse tridiagonal matrix (Julia v0.6) which works in matlab. When I try to factor L:

L = spdiagm((-1.*v,10.*v,-1.*v),(-1,0,1))
println(typeof(L))
cholfact(L)

I get:

SparseMatrixCSC{Float64,Int64}

Base.LinAlg.PosDefException(0)

Stacktrace:
 [1] #cholfact!#8(::Float64, ::Function, ::Base.SparseArrays.CHOLMOD.Factor{Float64}, ::Base.SparseArrays.CHOLMOD.Sparse{Float64}) at ./sparse/cholmod.jl:1360
 [2] (::Base.LinAlg.#kw##cholfact!)(::Array{Any,1}, ::Base.LinAlg.#cholfact!, ::Base.SparseArrays.CHOLMOD.Factor{Float64}, ::Base.SparseArrays.CHOLMOD.Sparse{Float64}) at ./<missing>:0
 [3] #cholfact#10(::Float64, ::Array{Int64,1}, ::Function, ::Base.SparseArrays.CHOLMOD.Sparse{Float64}) at ./sparse/cholmod.jl:1398
 [4] #cholfact#11(::Array{Any,1}, ::Function, ::SparseMatrixCSC{Float64,Int64}) at ./sparse/cholmod.jl:1438
 [5] cholfact(::SparseMatrixCSC{Float64,Int64}) at ./sparse/cholmod.jl:1438
 [6] include_string(::String, ::String) at ./loading.jl:515

However, the following works:

L2 = sparse(rand(5,5))
println(typeof(L2))
L3 = L'*L
println(typeof(L3))
cholfact(L3)

and both L2 and L3 are of type SparseMatrixCSC{Float64,Int64}.

2

I assume v = ones(n)? If yes, note that L[end,end] = 0 and thus your matrix is not SPD.

2 Likes
3

What is v? Please provide a minimum working example.

julia> v = fill(1.0, 3);
julia> L = spdiagm((-1.*v,10.*v,-1.*v),(-1,0,1));
julia> Ld = Array(L)
4×4 Array{Float64,2}:
 10.0  -1.0   0.0   0.0
 -1.0  10.0  -1.0   0.0
  0.0  -1.0  10.0  -1.0
  0.0   0.0  -1.0   0.0

Anyway, the problem is easier to see when it prints nicely (as a dense array).
The last diagonal element is zero; this matrix is not positive definite.

EDIT:
One suggestion:

julia> L = spdiagm((-1.*@view(v[2:end]),10.*v,-1.*@view(v[1:end-1])),(-1,0,1))
3×3 SparseMatrixCSC{Float64,Int64} with 7 stored entries:
  [1, 1]  =  10.0
  [2, 1]  =  -1.0
  [1, 2]  =  -1.0
  [2, 2]  =  10.0
  [3, 2]  =  -1.0
  [2, 3]  =  -1.0
  [3, 3]  =  10.0

julia> Array(L)
3×3 Array{Float64,2}:
 10.0  -1.0   0.0
 -1.0  10.0  -1.0
  0.0  -1.0  10.0
1 Like
4

Woah thank you… Sorry for the Q; can I delete this post since it’s not a Julia issue?

5

We don’t yet have an explicit policy (though are in the process of formulating one), but we generally try to avoid deleting posts unless absolutely necessary.