Add full(F::Factorization) for various kinds of factorization

base/linalg/bunchkaufman.jl

@@ -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

base/linalg/cholesky.jl

@@ -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)

base/linalg/factorization.jl

@@ -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

base/linalg/ldlt.jl

@@ -57,3 +57,12 @@ function A_ldiv_B!{T}(S::LDLt{T,SymTridiagonal{T}}, B::AbstractVecOrMat{T})
     end
     return B
 end
+
+## reconstruct the original matrix, which is tridiagonal
+function full(F::LDLt)
+    e = copy(F.data.ev)
+    d = copy(F.data.dv)
+    e .*= d[1:end-1]
+    d[2:end] += e .* F.data.ev
+    SymTridiagonal(d, e)
+end

base/linalg/lu.jl

@@ -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]),:]

doc/stdlib/linalg.rst

@@ -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``.

test/linalg1.jl

@@ -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")

test/linalg2.jl

@@ -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(full(Tldlt)) Fs
     end
 
     # eigenvalues/eigenvectors of symmetric tridiagonal