I think if you check *julia> **F.info*
*1*
Which is an error code (l.237 in LAPACK). I think this is because
*julia> **eigvals(A)*
*10-element Array{Float64,1}:*
* -2.97808e-16*
* -1.28759e-16*
* -5.1601e-17 *
* 2.89122e-17*
* 6.39045e-17*
* 1.47151e-16*
* 6.33471e-16*
* 0.184705 *
* 0.848542 *
* 8.02861 *
you have three negative eigenvalues (to machine precision). If you print L'
* L - P' * A * P, you can see where the factorization just stops.
If you perturb:
*julia> **B = A + eps() * eye(10)*
the same process gets you a difference very close to 0 everywhere (even
though it also exits with error code 1).
On Thursday, February 11, 2016 at 7:27:27 PM UTC-5, Michael Lindon wrote:
>
> The covariance matrix of my multivariate normal is not full rank and so
> regular cholesky throws a positive definite error.
>
> However, it is possible to do a pivoted cholesky decomposition. Quoting
> from the netlib manual of pstrf
> :
> http://www.netlib.org/lapack/explore-html/dd/dad/dpstrf_8f.html
>
> The factorization has the form
> P**T * A * P = U**T * U , if UPLO = 'U',
> P**T * A * P = L * L**T, if UPLO = 'L',
> where U is an upper triangular matrix and L is lower triangular, and
> P is stored as vector PIV.
>
>
>
> I want to generate a MvNormal draw by multiplying P*L*z where z is a
> vector of standard normals.
>
>
> From the lapack doc it seems I should be able to recover A from P*L*L'*P'.
> I have a minimial reproducible example where I try to recover A, but some
> of the elements of the difference matrix are larger than i would expect or
> hope:
>
> srand(1)
>
>
> a=rand(3,10)
>
>
> A=a'*a
>
>
>
>
>
> rank(A)
>
>
>
>
>
> F=cholfact(A,:L,Val{true})
>
>
> L=F[:L]
>
>
> P=F[:P]
>
>
> P*L*L'*P'-A
>
>
>
> julia> P*L*L'*P'-A
> 10x10 Array{Float64,2}:
> 1.07228 0.0 0.0 1.27149 1.11022e-16 0.0
> -5.55112e-17 1.84107 2.01638 0.0
> 0.0 0.0 0.0 0.0 0.0 0.0
> 0.0 0.0 0.0 0.0
> 0.0 0.0 0.0 -2.22045e-16 0.0 0.0
> 0.0 2.22045e-16 0.0 0.0
> 1.27149 0.0 -2.22045e-16 3.11518 -3.33067e-16 0.0
> 0.0 2.82245 2.48629 0.0
> 1.11022e-16 0.0 0.0 -3.33067e-16 5.55112e-16 0.0
> 0.0 0.0 2.22045e-16 0.0
> 0.0 0.0 0.0 0.0 0.0 0.0
> -5.55112e-17 0.0 0.0 0.0
> -5.55112e-17 0.0 0.0 0.0 0.0 -5.55112e-17
> 0.0 0.0 0.0 0.0
> 1.84107 0.0 2.22045e-16 2.82245 0.0 0.0
> 0.0 4.15206 3.59957 0.0
> 2.01638 0.0 0.0 2.48629 2.22045e-16 0.0
> 0.0 3.59957 4.09992 0.0
> 0.0 0.0 0.0 0.0 0.0 0.0
> 0.0 0.0 0.0 0.0
>
>
> notice some elements are large. Why is it not identically zero?
>
>
