Skip to content

Commit

Permalink
Add full(F::Factorization) for various kinds of factorization
Browse files Browse the repository at this point in the history
Include changes suggested by Andreas Noack
  • Loading branch information
davidavdav committed Dec 10, 2014

Verified

This commit was created on GitHub.com and signed with GitHub’s verified signature. The key has expired.
1 parent 54511c3 commit 1c11ed9
Showing 8 changed files with 37 additions and 0 deletions.
4 changes: 4 additions & 0 deletions base/linalg/bunchkaufman.jl
Original file line number Diff line number Diff line change
@@ -43,3 +43,7 @@ A_ldiv_B!{T<:BlasReal}(B::BunchKaufman{T}, R::StridedVecOrMat{T}) = LAPACK.sytrs
function A_ldiv_B!{T<:BlasComplex}(B::BunchKaufman{T}, R::StridedVecOrMat{T})
(issym(B) ? LAPACK.sytrs! : LAPACK.hetrs!)(B.uplo, B.LD, B.ipiv, R)
end

## reconstruct the original matrix
## TODO: understand the procedure described at
## http://www.nag.com/numeric/FL/nagdoc_fl22/pdf/F07/f07mdf.pdf
1 change: 1 addition & 0 deletions base/linalg/cholesky.jl
Original file line number Diff line number Diff line change
@@ -108,6 +108,7 @@ convert{T}(::Type{Factorization{T}}, C::CholeskyPivoted) = convert(CholeskyPivot

full{T,S}(C::Cholesky{T,S,:U}) = C[:U]'C[:U]
full{T,S}(C::Cholesky{T,S,:L}) = C[:L]*C[:L]'
full(F::CholeskyPivoted) = (ip=invperm(F[:p]); (F[:L] * F[:U])[ip,ip])

size(C::Union(Cholesky, CholeskyPivoted)) = size(C.UL)
size(C::Union(Cholesky, CholeskyPivoted), d::Integer) = size(C.UL,d)
14 changes: 14 additions & 0 deletions base/linalg/factorization.jl
Original file line number Diff line number Diff line change
@@ -716,3 +716,17 @@ function At_ldiv_B{TF<:Number,TB<:Number,N}(F::Factorization{TF}, B::AbstractArr
TFB = typeof(one(TF)/one(TB))
At_ldiv_B!(convert(Factorization{TFB}, F), TB == TFB ? copy(B) : convert(AbstractArray{TFB,N}, B))
end

## reconstruct the original matrix
full(F::QR) = F[:Q] * F[:R]
full(F::QRCompactWY) = F[:Q] * F[:R]
full(F::QRPivoted) = (F[:Q] * F[:R])[:,invperm(F[:p])]

## Can we determine the source/result is Real? This is not stored in the type Eigen
full(F::Eigen) = F.vectors * Diagonal(F.values) / F.vectors

full(F::Hessenberg) = (fq = full(F[:Q]); fq * F[:H] * fq')

full(F::Schur) = F.Z * F.T * F.Z'

full(F::SVD) = F.U * Diagonal(F.S) * F.Vt
3 changes: 3 additions & 0 deletions base/linalg/ldlt.jl
Original file line number Diff line number Diff line change
@@ -57,3 +57,6 @@ function A_ldiv_B!{T}(S::LDLt{T,SymTridiagonal{T}}, B::AbstractVecOrMat{T})
end
return B
end

## reconstruct the original matrix, but in full rather than tridiagonal
full(F::LDLt) = (L=Triangular(F.data, :L, true); L * Diagonal(diag(F.data)) * L')
2 changes: 2 additions & 0 deletions base/linalg/lu.jl
Original file line number Diff line number Diff line change
@@ -318,3 +318,5 @@ end

/(B::AbstractMatrix,A::LU) = At_ldiv_Bt(A,B).'

## reconstruct the original matrix
full(F::LU) = (F[:L] * F[:U])[invperm(F[:p]),:]
4 changes: 4 additions & 0 deletions doc/stdlib/linalg.rst
Original file line number Diff line number Diff line change
@@ -41,6 +41,10 @@ Linear algebra functions in Julia are largely implemented by calling functions f

``factorize!`` is the same as :func:`factorize`, but saves space by overwriting the input ``A``, instead of creating a copy.

.. function:: full(F)

Reconstruct the matrix ``A`` from the factorization ``F=factorize(A)``.

.. function:: lu(A) -> L, U, p

Compute the LU factorization of ``A``, such that ``A[p,:] = L*U``.
8 changes: 8 additions & 0 deletions test/linalg1.jl
Original file line number Diff line number Diff line change
@@ -78,6 +78,7 @@ debug && println("pivoted Choleksy decomposition")
if isreal(apd)
@test_approx_eq apd * inv(cpapd) eye(n)
end
@test_approx_eq full(cpapd) apd
end

debug && println("(Automatic) Bunch-Kaufman factor of indefinite matrix")
@@ -100,6 +101,7 @@ debug && println("(Automatic) Square LU decomposition")
@test_approx_eq (l*u)[invperm(p),:] a
@test_approx_eq a * inv(lua) eye(n)
@test norm(a*(lua\b) - b, 1) < ε*κ*n*2 # Two because the right hand side has two columns
@test_approx_eq full(lua) a

debug && println("Thin LU")
lua = lufact(a[:,1:5])
@@ -116,6 +118,7 @@ debug && println("QR decomposition (without pivoting)")
@test_approx_eq q*full(q, thin=false)' eye(n)
@test_approx_eq q*r a
@test_approx_eq_eps a*(qra\b) b 3000ε
@test_approx_eq full(qra) a

debug && println("(Automatic) Fat (pivoted) QR decomposition") # Pivoting is only implemented for BlasFloats
qrpa = factorize(a[1:5,:])
@@ -126,6 +129,7 @@ debug && println("(Automatic) Fat (pivoted) QR decomposition") # Pivoting is onl
@test_approx_eq q*r isa(qrpa,QRPivoted) ? a[1:5,p] : a[1:5,:]
@test_approx_eq isa(qrpa, QRPivoted) ? q*r[:,invperm(p)] : q*r a[1:5,:]
@test_approx_eq_eps a[1:5,:]*(qrpa\b[1:5]) b[1:5] 5000ε
@test_approx_eq full(qrpa) a[1:5,:]

debug && println("(Automatic) Thin (pivoted) QR decomposition") # Pivoting is only implemented for BlasFloats
qrpa = factorize(a[:,1:5])
@@ -135,6 +139,7 @@ debug && println("(Automatic) Thin (pivoted) QR decomposition") # Pivoting is on
@test_approx_eq q*full(q, thin=false)' eye(n)
@test_approx_eq q*r isa(qrpa, QRPivoted) ? a[:,p] : a[:,1:5]
@test_approx_eq isa(qrpa, QRPivoted) ? q*r[:,invperm(p)] : q*r a[:,1:5]
@test_approx_eq full(qrpa) a[:,1:5]

debug && println("symmetric eigen-decomposition")
if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia
@@ -146,6 +151,7 @@ debug && println("symmetric eigen-decomposition")
@test_approx_eq abs(eigfact(Hermitian(asym), d[1]-10*eps(d[1]), d[2]+10*eps(d[2]))[:vectors]'v[:,1:2]) eye(eltya, 2)
@test_approx_eq eigvals(Hermitian(asym), 1:2) d[1:2]
@test_approx_eq eigvals(Hermitian(asym), d[1]-10*eps(d[1]), d[2]+10*eps(d[2])) d[1:2]
@test_approx_eq full(eigfact(asym)) asym
end

debug && println("non-symmetric eigen decomposition")
@@ -178,6 +184,7 @@ debug && println("Schur")
@test_approx_eq sort(real(f[:values])) sort(real(d))
@test_approx_eq sort(imag(f[:values])) sort(imag(d))
@test istriu(f[:Schur]) || iseltype(a,Real)
@test_approx_eq full(f) a
end

debug && println("Reorder Schur")
@@ -205,6 +212,7 @@ debug && println("singular value decomposition")
if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia
usv = svdfact(a)
@test_approx_eq usv[:U]*scale(usv[:S],usv[:Vt]) a
@test_approx_eq full(usv) a
end

debug && println("Generalized svd")
1 change: 1 addition & 0 deletions test/linalg2.jl
Original file line number Diff line number Diff line change
@@ -82,6 +82,7 @@ for elty in (Float32, Float64, Complex64, Complex128, Int)
Tldlt = ldltfact(Ts)
x = Tldlt\v
@test_approx_eq x invFsv
@test_approx_eq full(Tldlt) Fs
end

# eigenvalues/eigenvectors of symmetric tridiagonal

0 comments on commit 1c11ed9

Please sign in to comment.