Added choleskyinv.jl unit

.travis.yml

@@ -5,8 +5,6 @@ os:
   - osx
   - windows
 julia:
-  - 1.0
-  - 1.1
   - 1.2
   - 1.3
   - nightly

Project.toml

@@ -1,6 +1,6 @@
 name = "MatrixFactorizations"
 uuid = "a3b82374-2e81-5b9e-98ce-41277c0e4c87"
-version = "0.2.1"
+version = "0.3.0"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

src/MatrixFactorizations.jl

@@ -8,7 +8,7 @@ import LinearAlgebra.LAPACK: liblapack, chkuplo, chktrans
 import LinearAlgebra: cholesky, cholesky!, norm, diag, eigvals!, eigvals, eigen!, eigen,
             qr, axpy!, ldiv!, mul!, lu, lu!, ldlt, ldlt!, AbstractTriangular,
             chkstride1, kron, lmul!, rmul!, factorize, StructuredMatrixStyle, logabsdet,
-            AbstractQ, _zeros, _cut_B, _ret_size
+            AbstractQ, _zeros, _cut_B, _ret_size, require_one_based_indexing, checksquare
 
 import Base: getindex, setindex!, *, +, -, ==, <, <=, >,
                 >=, /, ^, \, transpose, showerror, reindex, checkbounds, @propagate_inbounds
@@ -23,17 +23,13 @@ import Base: convert, size, view, unsafe_indices,
 
 import ArrayLayouts: reflector!, reflectorApply!
 
-if VERSION < v"1.2-"
-    import Base: has_offset_axes
-    require_one_based_indexing(A...) = !has_offset_axes(A...) || throw(ArgumentError("offset arrays are not supported but got an array with index other than 1"))
-else
-    import Base: require_one_based_indexing    
-end                            
+
 
-export ql, ql!, qrunblocked, qrunblocked!, QL                        
+export ql, ql!, qrunblocked, qrunblocked!, QL, choleskyinv!, choleskyinv                        
 
 
 include("qr.jl")
 include("ql.jl")
+include("choleskyinv.jl")
 
 end #module

src/choleskyinv.jl

@@ -0,0 +1,143 @@
+# choleskyinv.jl
+# Cholesky factor and its inverse (complex-conjugate) transposed in one pass.
+# For small matrices this is faster than computing the Cholesky factor
+# using the `cholesky` function in LinearAlgebra and then invert it.
+
+
+"""
+Cholesky factorization and inverse Cholesky factorization, accessible
+in fields `.c` ans `.ci`, respectively.
+"""
+struct CholeskyInv{T, F<:Factorization{T}} <: Factorization{T}
+     c::F
+     ci::F
+end
+
+"""
+    choleskyinv(P::AbstractMatrix{T};
+		check::Bool=true, tol::Real = √eps(T)) where T<:Union{Real, Complex}
+
+ Compute the Cholesky factorization of a dense positive definite
+ matrix P and return a `Choleskyinv` object, holding in field `.c`
+ the Cholesky factorization and in field `ci` the inverse of the Cholesky
+ factorization.
+
+ The two factorizations are obtained in one pass and this is faster
+ then calling Julia's [chelosky](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.cholesky)
+ function and inverting the lower factor for small matrices.
+
+ Input matrix `P` may be of type `Matrix` or `Hermitian`. Since only the
+ lower triangle is used, `P` may also be a `LowerTriangular` view of a
+ positive definite matrix.
+ If `P` is real, it can also be of the `Symmetric` type.
+
+ The algorithm is a *multiplicative Gaussian elimination*.
+
+ **Notes:**
+ The inverse Cholesky factor ``L^{-1}``, obtained as `.ci.L`,
+ is an inverse square root (whitening matrix) of `P`, since
+ ``L^{-1}PL^{-H}=I``. It therefore yields the inversion of ``P`` as
+ ``P^{-1}=L^{-H}L^{-1}``. This is the fastest whitening matrix to be computed,
+ however it yields poor numerical precision, especially for large matrices.
+
+ The following relations holds:
+ - ``L=PL^{-H}``,
+ - ``L^{H}=L^{-1}P``,
+ - ``L^{-H}=P^{-1}L``
+ - ``L^{-1}=L^{H}P^{-1}``.
+
+ We also have ``L^{H}L=L^{-1}P^{2}L^{-H}=UPU^H``, with ``U`` orthogonal
+ (see below) and ``L^{-1}L^{-H}=L^{H}P^{-2}L=UP^{-1}U^H``.
+ ``LL^{H}`` and ``L^{H}L`` are unitarily similar, that is,
+
+ ``ULL^{H}U^H=L^{H}L``,
+
+ where ``U=L^{-1}P^{1/2}``, with ``P^{1/2}=H`` the *principal* (unique symmetric)
+ square root of ``P``. This is seen writing
+ ``PP^{-1}=HHL^{-H}L^{-1}``; multiplying both sides on the left by ``L^{-1}``
+ and on the right by ``L`` we obtain
+ ``L^{-1}PP^{-1}L=L^{-1}HHL^{-H}=I=(L^{-1}H)(L^{-1}H)^H`` and since
+ ``L^{-1}H`` is square it must be unitary.
+
+ From these expressions we have
+ - ``H=LU=U^HL^H,
+ - ``L=HU^H``,
+ - ``H^{-1}=U^HL^{-1}``
+ - ``L^{-1}=UH^{-1}``.
+
+ ``U`` is the *polar factor* of ``L^{H}``, *i.e.*, ``L^{H}=UH``,
+ since ``LL^{H}=HU^HUH^H=H^2=P``.
+ From ``L^{H}L=UCU^H`` we have ``L^{H}LU=UC=ULL^{H}`` and from
+ ``U=L^{-1}H`` we have ``L=HU^H``.
+
+ ## Examples
+	using PosDefManifold, Test
+	n = 40
+	etol = 1e-9
+	Y=randP(n)
+	Yi=inv(Y)
+	a=cholesky(Y)
+
+	C=choleskyinv!(copy(Matrix(Y)))
+	@test(norm(C.c.L*C.c.U-Y)/√n < etol)
+	@test(norm(C.ci.U*C.ci.L-Yi)/√n < etol)
+
+	# repeat the test for complex matrices
+	Y=randP(ComplexF64, n)
+	Yi=inv(Y)
+	C=choleskyinv!(copy(Matrix(Y)))
+	@test(norm(C.c.L*C.c.U-Y)/√n < etol)
+	@test(norm(C.ci.U*C.ci.L-Yi)/√n < etol)
+
+	# Benchmark
+	using BenchmarkTools
+
+	# computing the Cholesky factor and its inverse using LinearAlgebra
+	function linearAlgebraWay(P)
+		C=cholesky(P)
+		Q=inv(C.L)
+	end
+
+	Y=randP(n)
+	@benchmark(choleskyinv(Y))
+	@benchmark(linearAlgebraWay(Y))
+
+"""
+choleskyinv(P::Union{Hermitian, Symmetric, Matrix, LowerTriangular};
+	   		check::Bool = true,
+	   		tol::Real = √eps(real(eltype(P)))) =
+    choleskyinv!(copy(Matrix(P)); check=check, tol=tol)
+
+"""
+    choleskyinv!(P::AbstractMatrix{T};
+		kind::Symbol = :LLt, tol::Real = √eps(T)) where T<:RealOrComplex
+ The same thing as [`choleskyinv`](@ref), but destroys the input matrix.
+"""
+function choleskyinv!(P::Matrix{T};
+			check::Bool = true,
+			tol::Real = √eps(real(T))) where T<:Union{Real, Complex}
+	LinearAlgebra.require_one_based_indexing(P)
+	n = LinearAlgebra.checksquare(P)
+
+	@inbounds for j=1:n-1
+		check && abs2(P[j, j])<tol && throw(LinearAlgebra.PosDefException(1))
+		for i=j+1:n
+			θ = conj(P[i, j] / -P[j, j])
+			for k=i:n P[k, i] += θ * P[k, j] end # update P and write D
+			P[i, j] = conj(-θ) # write Cholesky factor
+			for k=1:j-1 P[k, i] += θ * P[k, j] end # write inv Cholesky factor
+			P[j, i] = θ # write inv Cholesky factor
+		end
+	end
+
+	LT, ULT, UUT = LowerTriangular, UnitLowerTriangular, UnitUpperTriangular
+	D=sqrt.(Diagonal(P))
+	return CholeskyInv(Cholesky(LT(ULT(P)*D), :L, 0), Cholesky(LT(Matrix((UUT(P)*inv(D))')) , :L, 0 ))
+end
+
+
+function show(io::IO, ::MIME{Symbol("text/plain")}, C::CholeskyInv)
+    println(io, "Cholesky (.c) and inverse Cholesky (.ci) factorizations")
+	#show(io, C.c)
+	#show(io, C.ci)
+end

test/runtests.jl

@@ -1,5 +1,5 @@
 using MatrixFactorizations, LinearAlgebra, Random, Test
-using LinearAlgebra: BlasComplex, BlasFloat, BlasReal, rmul!, lmul!
+using LinearAlgebra: BlasComplex, BlasFloat, BlasReal, rmul!, lmul!, require_one_based_indexing, checksquare
 
 n = 10
 
@@ -378,4 +378,31 @@ end
         @test rmul!(copy(c), Q) ≈ c*Matrix(Q)
         @test rmul!(copy(c), Q') ≈ c*Matrix(Q')
     end
-end
+end
+
+# choleskyinv.jl
+@testset "choleskyinv" begin
+    etol = 1e-9
+    n = 20
+	# real matrices
+    A = randn(n, n)
+    P = A*A'
+    Pi = inv(P)
+    C = choleskyinv!(copy(Matrix(P)))
+	@test(norm(C.c.L*C.c.U-P)/√n < etol)
+	@test(norm(C.ci.U*C.ci.L-Pi)/√n < etol)
+    C = choleskyinv(P)
+	@test(norm(C.c.L*C.c.U-P)/√n < etol)
+	@test(norm(C.ci.U*C.ci.L-Pi)/√n < etol)
+
+	# repeat the test for complex matrices
+    A=randn(ComplexF64, n, n)
+    P=A*A'
+    Pi=inv(P)
+	C=choleskyinv!(copy(Matrix(P)))
+	@test(norm(C.c.L*C.c.U-P)/√n < etol)
+	@test(norm(C.ci.U*C.ci.L-Pi)/√n < etol)
+    C = choleskyinv(P)
+	@test(norm(C.c.L*C.c.U-P)/√n < etol)
+	@test(norm(C.ci.U*C.ci.L-Pi)/√n < etol)
+end