juliastats / pdmats.jl Goto Github PK

It shows error in lu.jl file at following line-
function lufact{T}(A::AbstractMatrix{T}, pivot::Union{Type{Val{false}}, Type{Val{true}}})
S = typeof(zero(T)/one(T))
AA = similar(A, S, size(A))
copy!(AA, A)
lufact!(AA, pivot)
end

It shows following error
LoadError: MethodError: no method matching zero(::Type{Any})
Closest candidates are:
zero(!Matched::Type{Base.LibGit2.Oid}) at libgit2\oid.jl:88
zero(!Matched::Type{Base.Pkg.Resolve.VersionWeights.VWPreBuildItem}) at pkg/resolve/versionweight.jl:80
zero(!Matched::Type{Base.Pkg.Resolve.VersionWeights.VWPreBuild}) at pkg/resolve/versionweight.jl:120
...
in lufact(::Array{Any,2}) at .\linalg\lu.jl:118
in det(::Array{Any,2}) at .\linalg\generic.jl:578
in disp_bwo(::Float64, ::Complex{Float64}) at C:\Users\admin\Desktop\new julia program\runit1.jl:208
in search_roots(::Float64, ::Array{Float64,1}) at C:\Users\admin\Desktop\new julia program\runit1.jl:316
in include_string(::String, ::String) at .\loading.jl:441
in include_string(::Module, ::String, ::String) at C:\Users\admin.julia\v0.5\CodeTools\src\eval.jl:32
in (::LastMain.LastMain.Atom.##61#64{String,String})() at C:\Users\admin.julia\v0.5\Atom\src\eval.jl:81
in withpath(::LastMain.LastMain.Atom.##61#64{String,String}, ::String) at C:\Users\admin.julia\v0.5\CodeTools\src\utils.jl:30
in withpath(::Function, ::String) at C:\Users\admin.julia\v0.5\Atom\src\eval.jl:46
in macro expansion at C:\Users\admin.julia\v0.5\Atom\src\eval.jl:79 [inlined]
in (::LastMain.LastMain.Atom.##60#63{String,String})() at .\task.jl:60
while loading C:\Users\admin\Desktop\new julia program\runit1.jl, in expression starting on line 389

julia> a = PDMat(eye(3))
PDMats.PDMat{Float64,Array{Float64,2}}(3,[1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0],Base.LinAlg.Cholesky{Float64,Array{Float64,2}} with factor:
[1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0])

julia> 3 * a
ERROR: StackOverflowError:
 in *(::PDMats.PDMat{Float64,Array{Float64,2}}, ::Int64) at /Users/cbinz/.julia/v0.5/PDMats/src/generics.jl:24 (repeats 80000 times)

This works as expected when the scaling factor is a Float64, by the way.

The UniformScaling type provides a convenient notation for a ScalMat in cases where the dimension can be inferred. By defining methods like

+(a::AbstractPDMat,b::UniformScaling) = a + ScalMat(dim(a), @compat(Float64(b.λ)))

a common operation like adding a multiple of the identity can be written more concisely.

I can provide a PR if this seems like a good idea.

In this line:

I ran into this issue doing the following:

using LinearAlgebra
using Distributions
using PDMats

J = sprand(10, 10, 0.05)
J = J'J
J[diagind(J)] .= 1
μ = randn(10)
h = J \ μ
d = MvNormalCanon(μ, h, PDSparseMat(J))
rand(d)

N.B., I also don't think sampling from an MvNormalCanon should require full inversion of the precision matrix, but that's a separate issue for Distributions.

Thanks for providing this extremely useful package!

I was wondering if it would be reasonable to include sparse covariance matrices. I'd be happy to help, although my knowledge about sparse algorithm is very limited.

quad(a, x)          # compute x' * a * x when `x` is a vector.
X_A_Xt(a, x)        # compute `x * a * x'` for a matrix `x`.

Right now, to take advantage of quad and X_A_Xt the user must ask for them directly.
Which is annoying because the user might not know they are dealing with a PDMat.
If the user wrote x' * a * x then it would not hit these methods.

But we could fix this via something like:

function Base.(:*)(x::Union{Adjoint{<:AbstractVector}, Transpose{<:AbstractVector},  a::AbstractPDMat, y::AbstractVector)
    if partent(x) == y
        quad(a, x)
    else
        invoke(...)  # call the generic `*`
    end
end

And similar for X_A_Xt

Hi!

It would be nice to have a new release for PDMats. Thank you.

julia> PDMat(eye(2)) * ones(2,1)
2x1 Array{Float64,2}:
 1.0
 1.0

julia> ones(1,2) * PDMat(eye(2))
ERROR: `*` has no method matching *(::Array{Float64,2}, ::PDMat)

(of course, it is possible to do (PDMat(eye(2))' * ones(1,2)' ) ')

Is this intentional?

Only the right division is tested right now

PackageEvaluator.jl is a script that runs nightly. It attempts to load all Julia packages and run their tests (if available) on both the stable version of Julia (0.3) and the nightly build of the unstable version (0.4). The results of this script are used to generate a package listing enhanced with testing results.

On Julia 0.4

• On 2014-10-27 the testing status was Package doesn't load.
• On 2014-10-28 the testing status changed to Tests fail, but package loads.

Tests fail, but package loads. means that PackageEvaluator found the tests for your package, executed them, and they didn't pass. However, trying to load your package with using worked.

Package doesn't load. means that PackageEvaluator did not find tests for your package. Additionally, trying to load your package with using failed.

This issue was filed because your testing status became worse. No additional issues will be filed if your package remains in this state, and no issue will be filed if it improves. If you'd like to opt-out of these status-change messages, reply to this message saying you'd like to and @IainNZ will add an exception. If you'd like to discuss PackageEvaluator.jl please file an issue at the repository. For example, your package may be untestable on the test machine due to a dependency - an exception can be added.

Test log:

>>> 'Pkg.add("PDMats")' log
INFO: Installing ArrayViews v0.4.8
INFO: Installing PDMats v0.2.4
INFO: Package database updated
INFO: METADATA is out-of-date a you may not have the latest version of PDMats
INFO: Use `Pkg.update()` to get the latest versions of your packages

>>> 'using PDMats' log
Julia Version 0.4.0-dev+1330
Commit 7fdc860 (2014-10-28 03:56 UTC)
Platform Info:
  System: Linux (x86_64-unknown-linux-gnu)
  CPU: Intel(R) Xeon(R) CPU E5-2650 0 @ 2.00GHz
  WORD_SIZE: 64
  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Sandybridge)
  LAPACK: libopenblas
  LIBM: libopenlibm
  LLVM: libLLVM-3.3

>>> test log
ERROR: type Cholesky has no field uplo
 in Xt_A_X at /home/idunning/pkgtest/.julia/v0.4/PDMats/src/PDMats.jl:222
 in include at ./boot.jl:242
 in include_from_node1 at ./loading.jl:128
 in include at ./boot.jl:242
 in include_from_node1 at loading.jl:128
 in process_options at ./client.jl:293
 in _start at ./client.jl:362
 in _start_3B_3769 at /home/idunning/julia04/usr/bin/../lib/julia/sys.so
while loading /home/idunning/pkgtest/.julia/v0.4/PDMats/test/pdmat.jl, in expression starting on line 98
while loading /home/idunning/pkgtest/.julia/v0.4/PDMats/test/runtests.jl, in expression starting on line 3
Running tests ...

INFO: Testing PDMats
===============================[ ERROR: PDMats ]================================

failed process: Process(`/home/idunning/julia04/usr/bin/julia /home/idunning/pkgtest/.julia/v0.4/PDMats/test/runtests.jl`, ProcessExited(1)) [1]

================================================================================
INFO: No packages to install, update or remove
ERROR: PDMats had test errors
 in error at error.jl:21
 in test at pkg/entry.jl:719
 in anonymous at pkg/dir.jl:28
 in cd at ./file.jl:20
 in cd at pkg/dir.jl:28
 in test at pkg.jl:68
 in process_options at ./client.jl:221
 in _start at ./client.jl:362
 in _start_3B_3769 at /home/idunning/julia04/usr/bin/../lib/julia/sys.so


>>> end of log

If we're not waiting on anything, could I request a new release?

I contributed #103 so that I could address some issues (mentioned here and here) with the matrix-variate distributions in Distributions.jl. Once these changes are released I can get that done.

For licensing issues we run a build of julia excluding GPL libraries. (USE_GPL_LIBS=0)

PDMats seems to be incompatible with this:

ERROR: LoadError: LoadError: UndefVarError: CHOLMOD not defined
 in include at ./boot.jl:261
 in include_from_node1 at ./loading.jl:304
 in include at ./boot.jl:261
 in include_from_node1 at ./loading.jl:304
while loading /dynizer/.julia_pkg/v0.4/PDMats/src/chol.jl, in expression starting on line 16
while loading /dynizer/.julia_pkg/v0.4/PDMats/src/PDMats.jl, in expression starting on line 50

The Cholesky bug is still breaking the travis tests on Distributions.jl. If you make a new release, then we can bump the requirements, which should fix it.

StridedMatrix is a annoying impossible extend union type.
it is not too hard for a custom array type to implement all the required methods.
(especially if it is a wrapper type that just delegates them all),

However it is not possible for a new type to join a Union.

I suggest that all dispatches on StridedMatrix in the package be replaced with AbstractMatrix.
If that throws an error deeper in the code because they hit some method that actually needs a StridedMatrix and they don't have the methods needed then so be it.

For A a PDMat and z a Vector, would it not be more numerically stable to compute invquad(A, z) as

zhalf  = A.chol.L \ z
q = dot(zhalf, zhalf)

rather than the current:

In some experiments my suggested way of computing invquad(A,z) more reliably gives positive values (as it always should!) when A and z contain very large values and A is close to singular. This is important when computing logpdf() of a multivariate normal distribution at fairly extreme values using the Distributions.jl package, which appears to use the PDMats invquad() function.

If I understand the source code correctly this is already how the computation is handled in the case of z being a matrix:

An analogous suggestion could also pertain to the computing of quad(A, z), though the numerical stability properties there are probably less of a problem. I'm running Julia 1.1.0 and PDMats 0.9.7.

Working with PDMats, it made me wonder why AbstractPDMat doesn't derive from AbstractMatrix? Given the class is immutable, we'd only have to implement getindex and size I think.

Is it possible to make a new release? I would like to make use of the bugfixes in #124.

LinearAlgebra.mul! returns an invalid result for PDiagMat. Given

using LinearAlgebra, PDMats

A = PDiagMat([0.25, 0.75])
x = [0.0, 1.0]
y = similar(x)

We get

A * x == [0.0, 0.75]

as expected, but mul! gives

mul!(y, A, x) == [0.0, 0.0]

However,

mul!(y, Matrix(A), x) == [0.0, 0.75]

returns the correct result.

For some reason, LinearAlgebra.generic_matvecmul! doesn't seem to work correctly for PDiagMat.

Tested with Julia version v1.5.3 on Linux x86_64.

CC @vindex10, who noticed incorrect results caused by this in a actual use case.

Ran into this while looking at autodiff for the MvNormal pdf, which calls sqmahal, which calls quad, which doesn't accept an array and vector with different eltypes.

I would propose replacing type signatures in utils.jl like

_addscal!{T<:Real}(r::Matrix{T}, a::Matrix{T}, b::Union{Matrix{T}, SparseMatrixCSC{T}}, c::T)

with

_addscal!(r::Matrix, a::Matrix, b::Union{Matrix, SparseMatrixCSC}, c::Real)

and making sure non-inplace versions in generics.jl allocate eltype(r) to be promote_type(eltype(a), eltype(b), eltype(c)) (perhaps with an @inplace function?). There's plenty of other code that simply duck-types the storage arrays, assuming their types are wide enough to hold the result of the computation, so this doesn't seem to be entirely out of place, but I wanted to make sure I wasn't missing something performance-critical here.

Happy to work on this PR, since I have to implement some version of this anyway.

PackageEvaluator.jl is a script that runs nightly. It attempts to load all Julia packages and run their tests (if available) on both the stable version of Julia (0.2) and the nightly build of the unstable version (0.3). The results of this script are used to generate a package listing enhanced with testing results.

On Julia 0.3

• On 2014-07-12 the testing status was Tests pass.
• On 2014-07-14 the testing status changed to Tests fail, but package loads.

Tests pass. means that PackageEvaluator found the tests for your package, executed them, and they all passed.

Tests fail, but package loads. means that PackageEvaluator found the tests for your package, executed them, and they didn't pass. However, trying to load your package with using worked.

This issue was filed because your testing status became worse. No additional issues will be filed if your package remains in this state, and no issue will be filed if it improves. If you'd like to opt-out of these status-change messages, reply to this message saying you'd like to and @IainNZ will add an exception. If you'd like to discuss PackageEvaluator.jl please file an issue at the repository. For example, your package may be untestable on the test machine due to a dependency - an exception can be added.

Test log:

INFO: Installing ArrayViews v0.4.6
INFO: Installing PDMats v0.2.0
INFO: Package database updated
Warning: could not import Base.add! into PDMats
Warning: could not import Base.add! into PDMats
ERROR: add! not defined
 in include at ./boot.jl:245
 in include_from_node1 at ./loading.jl:128
 in include at ./boot.jl:245
 in include_from_node1 at loading.jl:128
 in process_options at ./client.jl:285
 in _start at ./client.jl:354
while loading /home/idunning/pkgtest/.julia/v0.3/PDMats/test/pdmat.jl, in expression starting on line 248
while loading /home/idunning/pkgtest/.julia/v0.3/PDMats/runtests.jl, in expression starting on line 3
INFO: Package database updated

Note this is possibly due to removal of deprecated functions in Julia 0.3-rc1: JuliaLang/julia#7609

PackageEvaluator.jl is a script that runs nightly. It attempts to load all Julia packages and run their tests (if available) on both the stable version of Julia (0.3) and the nightly build of the unstable version (0.4). The results of this script are used to generate a package listing enhanced with testing results.

On Julia 0.4

• On 2014-08-14 the testing status was Tests pass.
• On 2014-08-15 the testing status changed to Tests fail, but package loads.

Tests pass. means that PackageEvaluator found the tests for your package, executed them, and they all passed.

Tests fail, but package loads. means that PackageEvaluator found the tests for your package, executed them, and they didn't pass. However, trying to load your package with using worked.

This issue was filed because your testing status became worse. No additional issues will be filed if your package remains in this state, and no issue will be filed if it improves. If you'd like to opt-out of these status-change messages, reply to this message saying you'd like to and @IainNZ will add an exception. If you'd like to discuss PackageEvaluator.jl please file an issue at the repository. For example, your package may be untestable on the test machine due to a dependency - an exception can be added.

Test log:

>>> 'Pkg.add("PDMats")' log
INFO: Installing ArrayViews v0.4.6
INFO: Installing PDMats v0.2.3
INFO: Package database updated

>>> 'using PDMats' log

>>> test log
Running tests ...

ERROR: type Cholesky has no field uplo
 in unwhiten at /home/idunning/pkgtest/.julia/v0.4/PDMats/src/PDMats.jl:180
 in include at ./boot.jl:245
 in include_from_node1 at ./loading.jl:128
 in include at ./boot.jl:245
 in include_from_node1 at loading.jl:128
 in process_options at ./client.jl:285
 in _start at ./client.jl:354
 in _start_3B_1724 at /home/idunning/julia04/usr/bin/../lib/julia/sys.so
while loading /home/idunning/pkgtest/.julia/v0.4/PDMats/test/pdmat.jl, in expression starting on line 59
while loading /home/idunning/pkgtest/.julia/v0.4/PDMats/test/runtests.jl, in expression starting on line 3

>>> end of log

Rather than linking to docs, the documentation link currently points to pkg.julialang.org. Even stranger, when I search there and eventually find PDMats, I get this:

I was wondering if it was remotely possible to use view on PDMat matrices. Technically one would just have to use view on the chol field and it should work out directly. However I am not sure what else would need to be implemented

At present there is no check for positivity of the arguments that will be on the diagonal of a PDiagMat or ScalMat. Probably a good idea to check that.

julia> using PDMats

julia> ScalMat(4,-1.5)
ScalMat(4,-1.5,-0.6666666666666666)

julia> PDiagMat([-1.,2.,-3.,4.])
PDiagMat(4,[-1.0,2.0,-3.0,4.0],[-1.0,0.5,-0.333333,0.25])

Currently, we get

julia> isposdef(PDMat([1.0 0.0; 0.0 1.0]))
ERROR: MethodError: no method matching cholesky(::Hermitian{Float64,PDMat{Float64,Array{Float64,2}}}; check=false)

julia> isposdef(PDiagMat([1.0, 1.0]))
ERROR: MethodError: no method matching cholesky(::Hermitian{Float64,PDiagMat{Float64,Array{Float64,1}}}; check=false)

julia> isposdef(ScalMat(2, 3.0))
ERROR: MethodError: no method matching cholesky(::Hermitian{Float64,ScalMat{Float64}}; check=false)

Shouldn't we define

LinearAlgebra.isposdef(::AbstractPDMat) = true

since AbstractPDMats are positive definite by definition?

Using version [90014a1f] PDMats v0.10.0 of PDMats when PDMats.isposdef called on a variable of type PDMat it raises a MethodError.
To reproduce

julia> s = [1. 0; 0 2]
2×2 Array{Float64,2}:
 1.0  0.0
 0.0  2.0

julia> PDMats.isposdef(s)
true

julia> PDMats.isposdef(PDMat(s))
ERROR: MethodError: no method matching cholesky(::LinearAlgebra.Hermitian{Float64,PDMat{Float64,Array{Float64,2}}}; check=false)

Possibly related to #117

PDMats is currently broken on Julia v0.7, since Base.Test was removed. This doesn't only affect PDMats testing but PDMats as a whole, as Base.Test is used in "src/testutils.jl".

It would be usedful if one could call cholesky on a AbstractPDMat

For PDMats this would just return the precomputed one,
For DiagPDMats this would just do the elementwise sqrt (right?).

cc @nickrobinson251

Dear all,

Is there a reason why I cannot do the following:

A = sparse([1,2],[1,2],[1.0,1.0])
B = Symmetric(A)
C= PDSparseMat(B)

This gives the following error:

MethodError: Cannot convert an object of type Symmetric{Float64,SparseMatrixCSC{Float64,Int64}} to an object of type PDMats.PDSparseMat

Best wishes,

Jan

Hi. First of all, this is a great idea and an important package to have for efficiency purposes.

Just a quick question/formality --I think it should also be specified that these matrices are hermitian(symmetric) in order for the Cholesky decomposition to be meaningful (see wiki). I could just be unfamiliar with the uses of the Cholesky decomposition for positive definite nonsymmetric matrices, but wanted to mention it just in case.

Also, the majority of your cases seem like they will be covariance matrices which are symmetric so this isn't relevant in those cases, but could possibly matter in cases like this where I don't think that scale matrices for an Inverse Wishart need be symmetric.

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

Is this just trivial to do?

@scheidan: Would you please look into this?

Testing PDMats.PDSparseMat with dim = 3
ERROR: LoadError: LoadError: assertion failed: |X_A_Xt(C,Xt) - Xt * Cmat * X| <= 4.440892098500626e-10
  X_A_Xt(C,Xt) = [1.316968930347459 0.18945306914517301 -1.4858054118599842 1.6380952858093882 1.249135489578486
 0.18945306914517301 6.4637003842729595 -3.458785536128513 -1.609462105240331 7.744455596263222
 -1.4858054118599842 -3.458785536128513 3.8318134295156163 -1.8071977928328986 -5.062616653189969
 1.6380952858093882 -1.609462105240331 -1.8071977928328986 4.088987807905241 -0.8897179122762191
 1.249135489578486 7.744455596263222 -5.062616653189969 -0.8897179122762191 10.125292113456553]
  Xt * Cmat * X = [0.2971437538066637 -0.053342197359836585 -0.16829829872370428 0.19428262828457893 0.20553115133824565
 -0.053342197359836585 1.0920671389496301 -0.5982762375146454 -0.20355504703958216 1.209785744990079
 -0.16829829872370428 -0.5982762375146454 0.6251578462929138 -0.2944772120091399 -0.7892472977327798
 0.19428262828457893 -0.20355504703958216 -0.2944772120091399 0.6367266496684335 -0.14715562747675964
 0.20553115133824562 1.209785744990079 -0.7892472977327798 -0.1471556274767596 1.5936973339384122]
  difference = 8.53159477951814 > 4.440892098500626e-10
 in error at error.jl:22
 in test_approx_eq at test.jl:140
 in test_approx_eq at test.jl:150
 in pdtest_triprod at /Users/dhlin/.julia/v0.4/PDMats/src/testutils.jl:201
 in test_pdmat at /Users/dhlin/.julia/v0.4/PDMats/src/testutils.jl:45
 in anonymous at no file:15
 in include at ./boot.jl:254
 in include_from_node1 at ./loading.jl:197
 in anonymous at no file:9
 in include at ./boot.jl:254
 in include_from_node1 at ./loading.jl:197
 in process_options at ./client.jl:308
 in _start at ./client.jl:411
while loading /Users/dhlin/.julia/v0.4/PDMats/test/pdmtypes.jl, in expression starting on line 8
while loading /Users/dhlin/.julia/v0.4/PDMats/test/runtests.jl, in expression starting on line 7
================================================================================[ ERROR: PDMats ]================================================================================

failed process: Process(`/Users/dhlin/julia/usr/bin/julia --check-bounds=yes --code-coverage=none --color=yes /Users/dhlin/.julia/v0.4/PDMats/test/runtests.jl`, ProcessExited(1)) [1]

=================================================================================================================================================================================
INFO: No packages to install, update or remove
ERROR: PDMats had test errors
 in error at ./error.jl:21
 in test at pkg/entry.jl:746
 in anonymous at pkg/dir.jl:31
 in cd at file.jl:22
 in cd at pkg/dir.jl:31
 in test at pkg.jl:71

whiten! and unwhiten use chol_lower to extract the lower triangular component of Cholesky decompositions. For SuiteSparse.CHOLMOD.Factor{T}, this extracts the lower triangular factor, but does not account for pivoting. This results in the same problem seen on this Discourse thread.

Using the sparse precision matrix attached to that thread,

using DataFrames, CSV

q_df = CSV.read("precision_matrix.csv")
point_df = CSV.read("mesh_points.csv")
Q = sparse(q_df.i, q_df.j, q_df.q)

Generate a realization,

using Random, SparseArrays, LinearAlgebra
using Distributions
using PDMats

Random.seed!(2)
z = randn(size(Q, 1))

pdQ = PDSparseMat(cholesky(Q))
xldiv = pdQ.chol.UP \ z
xwh = whiten(pdQ, z)

These two realizations should be the same, but they are not.

julia> maximum(xldiv - xwh)
45.00349329625478

This is also clear when they are plotted.

using Plots
pldiv = surface(point_df.x, point_df.y, xqs, camera=(45, 80), title="ldiv")
pwh = surface(point_df.x, point_df.y, xwh, camera=(45, 80), title="whiten")
plot(pldiv, pwh, layout=(1,2), size=(1200, 400), colorbar=false)

The ldiv version results in a random field with the expected structure, while the whiten version does not.

Downstream, this means that at least rand(::Distributions.MvNormalCanon) is giving incorrect results. I haven't looked closely at other methods in Distributions.jl yet.

Following up on the discussion in #34 it might be worth using the upper triangular. I shouldn't matter which one is used, but for type stability reasons we don't have an option for the lower version anymore on master.

Hi,

It seems that the invquad function is significantly slower on the master branch compared to the latest release v0.9.7.

using PDMats, LinearAlgebra, ForwardDiff, BenchmarkTools

A = PDMat(Matrix{Float64}(I, 250, 250));
x = rand(250);

@btime invquad($A, $x);

# v0.9.7
## 42.301 μs (1 allocation: 2.13 KiB)

# master
## 80.592 μs (3 allocations: 490.52 KiB)

### Duals ###

xd = ForwardDiff.Dual{Nothing}.(x)
@btime invquad($A, $xd);

# v0.9.7
## 40.114 μs (4 allocations: 2.17 KiB)

# master
## 330.757 μs (7 allocations: 978.94 KiB)

This slowdown is because the implementation of invquad for vectors went from:

invquad(M::PDMat, x) = dot(x, M.chol \ x)

to

invquad(M::PDMat, x) = sum(abs2, M.chol.L \ x)

According to the readme, this should work, and there is a method for it, but it fails

julia> c = Cholesky(UpperTriangular(q.R), :U, 0)
Cholesky{Float64,UpperTriangular{Float64,Array{Float64,2}}}
U factor:
4×4 UpperTriangular{Float64,Array{Float64,2}}:
 10.0084  -2.35472  -0.150555  -0.651091
   ⋅       8.48846  -0.832921   0.32629
   ⋅        ⋅       -9.36361    0.0278495
   ⋅        ⋅         ⋅        10.4634

julia> PDMat(c)
ERROR: MethodError: no method matching PDMat{Float64,Array{Float64,2}}(::Int64, ::Array{Float64,2}, ::Cholesky{Float64,UpperTriangular{Float64,Array{Float64,2}}})

[90014a1f] PDMats v0.9.11

The problem seems to be that the matrix in my cholesky object is an UpperTriangular as opposed to a regular matrix.

I installed a fresh julia 0.5.0 from FreeBSD ports. The default configuration compiles julia with USE_GPL_LIBS=0.

I then remove ~/.julia to start fresh:

[azo@euler ~]$ rm -rf .julia
[azo@euler ~]$ julia
               _
   _       _ _(_)_     |  A fresh approach to technical computing
  (_)     | (_) (_)    |  Documentation: http://docs.julialang.org
   _ _   _| |_  __ _   |  Type "?help" for help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 0.5.0 (2016-09-19 18:14 UTC)
 _/ |\__'_|_|_|\__'_|  |
|__/                   |  amd64-unknown-freebsd11.0

julia> Pkg.status()
INFO: Initializing package repository /home/azo/.julia/v0.5
INFO: Cloning METADATA from https://github.com/JuliaLang/METADATA.jl
No packages installed

julia> versioninfo()
Julia Version 0.5.0
Commit 3c9d753 (2016-09-19 18:14 UTC)
Platform Info:
  System: FreeBSD (amd64-unknown-freebsd11.0)
  CPU: Intel(R) Core(TM) i5-6500 CPU @ 3.20GHz
  WORD_SIZE: 64
  BLAS: libopenblas (NO_AFFINITY NEHALEM)
  LAPACK: libopenblas
  LIBM: libm
  LLVM: libLLVM-3.8.1 (ORCJIT, skylake)

julia>

At this point I am able to install PDMats but cannot do using PDMats:

julia> Pkg.add("PDMats")
INFO: Cloning cache of Compat from https://github.com/JuliaLang/Compat.jl.git
INFO: Cloning cache of PDMats from https://github.com/JuliaStats/PDMats.jl.git
INFO: Installing Compat v0.13.0
INFO: Installing PDMats v0.5.4
INFO: Package database updated

julia> using PDMats
INFO: Precompiling module PDMats.
ERROR: LoadError: LoadError: UndefVarError: CHOLMOD not defined
 in _init at /usr/local/lib/julia/sys.so:? (repeats 4 times)
 in macro expansion; at ./none:2 [inlined]
 in anonymous at ./<missing>:?
 in _init at /usr/local/lib/julia/sys.so:? (repeats 5 times)
while loading /home/azo/.julia/v0.5/PDMats/src/chol.jl, in expression starting on line 7
while loading /home/azo/.julia/v0.5/PDMats/src/PDMats.jl, in expression starting on line 50
ERROR: Failed to precompile PDMats to /home/azo/.julia/lib/v0.5/PDMats.ji.
 in _init at /usr/local/lib/julia/sys.so:? (repeats 5 times)
 in eval_user_input(::Any, ::Base.REPL.REPLBackend) at ./REPL.jl:64
 in macro expansion at ./REPL.jl:95 [inlined]
 in (::Base.REPL.##3#4{Base.REPL.REPLBackend})() at ./event.jl:68

From looking at chol.jl it seems that with julia v0.5.0 we should be using:

typealias CholTypeSparse{T} SparseArrays.CHOLMOD.Factor{T}

... but having compiled julia without USE_GPL_LIBS=1 I haven't got CHOLMOD in SparseArrays.
I tried using Base.LinAlg.CHOLMOD.CholmodFactor{T} instead but I realised I also don't have Base.LinAlg.CHOLMOD.

Therefore I decided to recompile julia with USE_GPL_LIBS=1. This compiled SparseArrays and included the necessary julia code and I could successfully install PDMats and do using PDMats

However, this time there were warnings when I started julia:

[azo@euler ~]$ rm -rf .julia
[azo@euler ~]$ julia
WARNING: Error during initialization of module CHOLMOD:
ErrorException("could not load library "libcholmod"
/usr/local/lib/libcholmod.so.1: Undefined symbol "colamd_printf"")
               _
   _       _ _(_)_     |  A fresh approach to technical computing
  (_)     | (_) (_)    |  Documentation: http://docs.julialang.org
   _ _   _| |_  __ _   |  Type "?help" for help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 0.5.0 (2016-09-19 18:14 UTC)
 _/ |\__'_|_|_|\__'_|  |
|__/                   |  amd64-unknown-freebsd11.0

julia> Pkg.add("PDMats")
INFO: Initializing package repository /home/azo/.julia/v0.5
INFO: Cloning METADATA from https://github.com/JuliaLang/METADATA.jl
INFO: Cloning cache of Compat from https://github.com/JuliaLang/Compat.jl.git
INFO: Cloning cache of PDMats from https://github.com/JuliaStats/PDMats.jl.git
INFO: Installing Compat v0.13.0
INFO: Installing PDMats v0.5.4
WARNING: Error during initialization of module CHOLMOD:
ErrorException("could not load library "libcholmod"
/usr/local/lib/libcholmod.so.1: Undefined symbol "colamd_printf"")
INFO: Package database updated

julia> using PDMats
INFO: Precompiling module PDMats.
WARNING: Error during initialization of module CHOLMOD:
ErrorException("could not load library "libcholmod"
/usr/local/lib/libcholmod.so.1: Undefined symbol "colamd_printf"")
WARNING: Error during initialization of module CHOLMOD:
ErrorException("could not load library "libcholmod"
/usr/local/lib/libcholmod.so.1: Undefined symbol "colamd_printf"")

julia>

So in essence I think there are two issues:

• isn't PDMats supposed to work with and without cholmod? if yes, why doesn't it work on FreeBSD?
• even with cholmod (USE_GPL_LIBS=1), we are we getting the "could not load library libcholmod" error?

I'm afraid I don't have a solution. /usr/local/lib/libcholmod.so.1 is available on my system but it doesn't seem to have colamd_printf in it.

Can someone please look into this?

I'm experimenting with using PDMats to represent covariance matrices for random-effects in the MixedModels package. Because I tend to use the factor more than the positive-definite matrix itself I will experiment with a different representation. When I am finished experimenting I will see how much overlap there is with this implementation and discuss whether to merge capabilities.

Base.LAPACK has been removed in Julia 0.7. However, Base.LinAlg.LAPACK still exists, even if it's deprecated.

The precompilation currently fails as Base.LAPACK cannot be imported on line 40 of PDMats.jl:

I would recommend using Base.LinAlg.LAPACK instead, which works for Julia < 0.7 and 0.7.

Ref: https://discourse.julialang.org/t/package-compatibility-caps/15301

PackageEvaluator.jl is a script that runs nightly. It attempts to load all Julia packages and run their tests (if available) on both the stable version of Julia (0.2) and the nightly build of the unstable version (0.3). The results of this script are used to generate a package listing enhanced with testing results.

On Julia 0.3

• On 2014-07-14 the testing status was Tests fail, but package loads.
• On 2014-07-15 the testing status changed to Package doesn't load.

Tests fail, but package loads. means that PackageEvaluator found the tests for your package, executed them, and they didn't pass. However, trying to load your package with using worked.

Package doesn't load. means that PackageEvaluator did not find tests for your package. Additionally, trying to load your package with using failed.

This issue was filed because your testing status became worse. No additional issues will be filed if your package remains in this state, and no issue will be filed if it improves. If you'd like to opt-out of these status-change messages, reply to this message saying you'd like to and @IainNZ will add an exception. If you'd like to discuss PackageEvaluator.jl please file an issue at the repository. For example, your package may be untestable on the test machine due to a dependency - an exception can be added.

Test log:

INFO: Installing ArrayViews v0.4.6
INFO: Installing PDMats v0.2.1
INFO: Package database updated
ERROR: NumericExtensions not found
 in require at loading.jl:47
 in include at ./boot.jl:245
 in include_from_node1 at ./loading.jl:128
 in reload_path at loading.jl:152
 in _require at loading.jl:67
 in require at loading.jl:51
 in include at ./boot.jl:245
 in include_from_node1 at loading.jl:128
 in process_options at ./client.jl:285
 in _start at ./client.jl:354
while loading /home/idunning/pkgtest/.julia/v0.3/PDMats/src/PDMats.jl, in expression starting on line 9
while loading /home/idunning/pkgtest/.julia/v0.3/PDMats/testusing.jl, in expression starting on line 1
INFO: Package database updated

Sorry for the elementary question! What I get is:

A = randn(5,5)
A = PDMat(A'*A)
B = randn(5,5)
B = PDMat(B'*B)
A*B
ERROR: MethodError: no method matching *(::PDMats.PDMat{Float64,Array{Float64,2}}, ::PDMats.PDMat{Float64,Array{Float64,2}})

Thanks in advance.

Shouldn't testutils.jl belong in the test folder instead of src ?

This is the sister issue of JuliaStats/Distributions.jl/issues/366 about normal distributions with degenerated covariance matrix. On implementation would be a famility of AbstractPSDMat's . The question is if this should also contain the concrete types now in the AbstractPDMat hierarchy. There are two natural choices of matrix factors, Cholesky with pivoting and LDL' decomposition, available form base respective from LAPACK.

Is there a specific reason why PDMats is not parameterised with respect to its floating point type?

Would it make sense to replace all instances of Float64 with T<:AbstractFloat, so that it can be used with Float32 instead of Float64?

The Cholesky will have eltype Float64 different from the than the original matrix eltype Int64, and the constructor does not allow that

julia> PDMat([1 0; 0 1])
ERROR: MethodError: no method matching PDMat{Int64,Array{Int64,2}}(::Int64, ::Array{Int64,2}, ::Cholesky{Float64,Array{Float64,2}})
Closest candidates are:
  PDMat{Int64,Array{Int64,2}}(::Int64, ::AbstractArray{T,2}, ::Cholesky{T,S}) where {T, S} at /Users/nick/.julia/packages/PDMats/AObTs/src/pdmat.jl:9
Stacktrace:
 [1] PDMat(::Array{Int64,2}, ::Cholesky{Float64,Array{Float64,2}}) at /Users/nick/.julia/packages/PDMats/AObTs/src/pdmat.jl:16
 [2] PDMat(::Array{Int64,2}) at /Users/nick/.julia/packages/PDMats/AObTs/src/pdmat.jl:19
 [3] top-level scope at none:0

Came up at JuliaStats/Distributions.jl#919 (comment)