Add Cholesky up- and downdates

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/base.rst

@@ -861,8 +861,8 @@ System
    * ``detach::Bool``\ : If ``true`` (defaults to ``false``\ ), then the ``Cmd`` will be   run in a new process group, allowing it to outlive the ``julia`` process   and not have Ctrl-C passed to it.
    * ``windows_verbatim::Bool``\ : If ``true`` (defaults to ``false``\ ), then on Windows   the ``Cmd`` will send a command-line string to the process with no quoting   or escaping of arguments, even arguments containing spaces.  (On Windows,   arguments are sent to a program as a single "command-line" string, and   programs are responsible for parsing it into arguments.  By default,   empty arguments and arguments with spaces or tabs are quoted with double   quotes ``"`` in the command line, and ``\`` or ``"`` are preceded by backslashes.   ``windows_verbatim=true`` is useful for launching programs that parse their   command line in nonstandard ways.)  Has no effect on non-Windows systems.
    * ``windows_hide::Bool``\ : If ``true`` (defaults to ``false``\ ), then on Windows no   new console window is displayed when the ``Cmd`` is executed.  This has   no effect if a console is already open or on non-Windows systems.
-   * ``env``\ : Set environment variables to use when running the ``Cmd``\ .  ``env``          is either a dictionary mapping strings to strings, an array          of strings of the form ``"var=val"``\ , an array or tuple of ``"var"=>val``          pairs, or ``nothing``\ .  In order to modify (rather than replace)          the existing environment, create ``env`` by ``copy(ENV)`` and then          set ``env["var"]=val`` as desired.
-   * ``dir::AbstractString``\ : Specify a working directory for the command (instead    of the current directory).
+   * ``env``\ : Set environment variables to use when running the ``Cmd``\ .  ``env``   is either a dictionary mapping strings to strings, an array   of strings of the form ``"var=val"``\ , an array or tuple of ``"var"=>val``   pairs, or ``nothing``\ .  In order to modify (rather than replace)   the existing environment, create ``env`` by ``copy(ENV)`` and then   set ``env["var"]=val`` as desired.
+   * ``dir::AbstractString``\ : Specify a working directory for the command (instead   of the current directory).
 
    For any keywords that are not specified, the current settings from ``cmd`` are used.   Normally, to create a ``Cmd`` object in the first place, one uses backticks, e.g.
 

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

doc/stdlib/pkg.rst

@@ -109,7 +109,7 @@ Functions for package development (e.g. ``tag``, ``publish``, etc.) have been mo
 
    .. Docstring generated from Julia source
 
-   Checkout the ``Pkg.dir(pkg)`` repo to the branch ``branch``\ . Defaults to checking out the "master" branch. To go back to using the newest compatible released version, use ``Pkg.free(pkg)``\ . Changes are merged (fast-forward only) if the keyword argument ``merge == true``, and the latest version is pulled from the upsream repo if ``pull == true``\ .
+   Checkout the ``Pkg.dir(pkg)`` repo to the branch ``branch``\ . Defaults to checking out the "master" branch. To go back to using the newest compatible released version, use ``Pkg.free(pkg)``\ . Changes are merged (fast-forward only) if the keyword argument ``merge == true``\ , and the latest version is pulled from the upsream repo if ``pull == true``\ .
 
 .. function:: pin(pkg)
 

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