Add Cholesky up- and downdates

Author: andreasnoack

base/linalg/cholesky.jl

@@ -264,3 +264,110 @@ end
 chkfullrank(C::CholeskyPivoted) = C.rank < size(C.factors, 1) && throw(RankDeficientException(C.info))
 
 rank(C::CholeskyPivoted) = C.rank
+
+"""
+    update!(C::Cholesky, v::StridedVector) -> CC::Cholesky
+
+Update a Cholesky factorization `C` with the vector `v`. If `A = C[:U]'C[:U]` then `CC = cholfact(C[:U]'C[:U] _ v*v')` but the computation of `CC` only uses `O(n^2)` operations. The input factorization `C` is updated in place such that on exit `C == CC`. The vector `v` is destroyed during the computation.
+"""
+function update!(C::Cholesky, v::StridedVector)
+    A = C.factors
+    n = length(v)
+    if size(C, 1) != n
+        throw(DimensionMismatch("updating vector must fit size of factorization"))
+    end
+    if C.uplo == 'U'
+        conj!(v)
+    end
+
+    for i = 1:n
+
+        # Compute Givens rotation
+        c, s, r = givensAlgorithm(A[i,i], v[i])
+
+        # Store new diagonal element
+        A[i,i] = r
+
+        # Update remaining elements in row/column
+        if C.uplo == 'U'
+            for j = i + 1:n
+                Aij = A[i,j]
+                vj  = v[j]
+                A[i,j]  =   c*Aij + s*vj
+                v[j]    = -s'*Aij + c*vj
+            end
+        else
+            for j = i + 1:n
+                Aji = A[j,i]
+                vj  = v[j]
+                A[j,i]  =   c*Aji + s*vj
+                v[j]    = -s'*Aji + c*vj
+            end
+        end
+    end
+    return C
+end
+
+"""
+    downdate!(C::Cholesky, v::StridedVector) -> CC::Cholesky
+
+Downdate a Cholesky factorization `C` with the vector `v`. If `A = C[:U]'C[:U]` then `CC = cholfact(C[:U]'C[:U] - v*v')` but the computation of `CC` only uses `O(n^2)` operations. The input factorization `C` is updated in place such that on exit `C == CC`. The vector `v` is destroyed during the computation.
+"""
+function downdate!(C::Cholesky, v::StridedVector)
+    A = C.factors
+    n = length(v)
+    if size(C, 1) != n
+        throw(DimensionMismatch("updating vector must fit size of factorization"))
+    end
+    if C.uplo == 'U'
+        conj!(v)
+    end
+
+    for i = 1:n
+
+        Aii = A[i,i]
+
+        # Compute Givens rotation
+        s = conj(v[i]/Aii)
+        s2 = abs2(s)
+        if s2 > 1
+            throw(LinAlg.PosDefException(i))
+        end
+        c = sqrt(1 - abs2(s))
+
+        # Store new diagonal element
+        A[i,i] = c*Aii
+
+        # Update remaining elements in row/column
+        if C.uplo == 'U'
+            for j = i + 1:n
+                vj = v[j]
+                Aij = (A[i,j] - s*vj)/c
+                A[i,j] = Aij
+                v[j] = -s'*Aij + c*vj
+            end
+        else
+            for j = i + 1:n
+                vj = v[j]
+                Aji = (A[j,i] - s*vj)/c
+                A[j,i] = Aji
+                v[j] = -s'*Aji + c*vj
+            end
+        end
+    end
+    return C
+end
+
+"""
+    update(C::Cholesky, v::StridedVector) -> CC::Cholesky
+
+Update a Cholesky factorization `C` with the vector `v`. If `A = C[:U]'C[:U]` then `CC = cholfact(C[:U]'C[:U] + v*v')` but the computation of `CC` only uses `O(n^2)` operations.
+"""
+update(C::Cholesky, v::StridedVector) = update!(copy(C), copy(v))
+
+"""
+    downdate(C::Cholesky, v::StridedVector) -> CC::Cholesky
+
+Downdate a Cholesky factorization `C` with the vector `v`. If `A = C[:U]'C[:U]` then `CC = cholfact(C[:U]'C[:U] - v*v')` but the computation of `CC` only uses `O(n^2)` operations.
+"""
+downdate(C::Cholesky, v::StridedVector) = downdate!(copy(C), copy(v))

doc/stdlib/linalg.rst

@@ -145,6 +145,34 @@ Linear algebra functions in Julia are largely implemented by calling functions f
 
    ``cholfact!`` is the same as :func:`cholfact`, but saves space by overwriting the input ``A``, instead of creating a copy. ``cholfact!`` can also reuse the symbolic factorization from a different matrix ``F`` with the same structure when used as: ``cholfact!(F::CholmodFactor, A)``.
 
+.. currentmodule:: Base.LinAlg
+
+.. function:: update(C::Cholesky, v::StridedVector) -> CC::Cholesky
+
+   .. Docstring generated from Julia source
+
+   Update a Cholesky factorization ``C`` with the vector ``v``\ . If ``A = C[:U]'C[:U]`` then ``CC = cholfact(C[:U]'C[:U] + v*v')`` but the computation of ``CC`` only uses ``O(n^2)`` operations.
+
+.. function:: downdate(C::Cholesky, v::StridedVector) -> CC::Cholesky
+
+   .. Docstring generated from Julia source
+
+   Downdate a Cholesky factorization ``C`` with the vector ``v``\ . If ``A = C[:U]'C[:U]`` then ``CC = cholfact(C[:U]'C[:U] - v*v')`` but the computation of ``CC`` only uses ``O(n^2)`` operations.
+
+.. function:: update!(C::Cholesky, v::StridedVector) -> CC::Cholesky
+
+   .. Docstring generated from Julia source
+
+   Update a Cholesky factorization ``C`` with the vector ``v``\ . If ``A = C[:U]'C[:U]`` then ``CC = cholfact(C[:U]'C[:U] _ v*v')`` but the computation of ``CC`` only uses ``O(n^2)`` operations. The input factorization ``C`` is updated in place such that on exit ``C == CC``\ . The vector ``v`` is destroyed during the computation.
+
+.. function:: downdate!(C::Cholesky, v::StridedVector) -> CC::Cholesky
+
+   .. Docstring generated from Julia source
+
+   Downdate a Cholesky factorization ``C`` with the vector ``v``\ . If ``A = C[:U]'C[:U]`` then ``CC = cholfact(C[:U]'C[:U] - v*v')`` but the computation of ``CC`` only uses ``O(n^2)`` operations. The input factorization ``C`` is updated in place such that on exit ``C == CC``\ . The vector ``v`` is destroyed during the computation.
+
+.. currentmodule:: Base
+
 .. function:: ldltfact(::SymTridiagonal) -> LDLt
 
    .. Docstring generated from Julia source

test/linalg/cholesky.jl

@@ -146,3 +146,13 @@ for elty in (Float32, Float64, Complex{Float32}, Complex{Float64})
     @test_approx_eq full(cholfact(A)[:L]) full(invoke(Base.LinAlg.chol!, Tuple{AbstractMatrix, Type{LowerTriangular}}, copy(A), LowerTriangular))
     @test_approx_eq full(cholfact(A)[:U]) full(invoke(Base.LinAlg.chol!, Tuple{AbstractMatrix, Type{UpperTriangular}}, copy(A), UpperTriangular))
 end
+
+# Test up- and downdates
+let A = complex(randn(10,5), randn(10, 5)), v = complex(randn(5), randn(5))
+    for uplo in (:U, :L)
+        AcA = A'A
+        F = cholfact(AcA, uplo)
+        @test LinAlg.update(F, v)[uplo] ≈ cholfact(AcA + v*v')[uplo]
+        @test LinAlg.downdate(F, v)[uplo] ≈ cholfact(AcA - v*v')[uplo]
+    end
+end