Add Cholesky up- and downdates

Author: andreasnoack

base/linalg.jl

@@ -126,6 +126,8 @@ export
     triu,
     tril!,
     triu!,
+    update,
+    update!,
     vecdot,
     vecnorm,
 

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

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