From 0b25019ca1ba043c13a8bf16924877eef3adee9c Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Tue, 1 Dec 2015 21:13:40 -0500 Subject: [PATCH] Add Cholesky up- and downdates --- base/linalg/cholesky.jl | 107 ++++++++++++++++++++++++++++++++++++++++ doc/stdlib/base.rst | 4 +- doc/stdlib/linalg.rst | 28 +++++++++++ doc/stdlib/pkg.rst | 2 +- test/linalg/cholesky.jl | 10 ++++ 5 files changed, 148 insertions(+), 3 deletions(-) diff --git a/base/linalg/cholesky.jl b/base/linalg/cholesky.jl index b75429d8723e9b..3604cdce9aa612 100644 --- a/base/linalg/cholesky.jl +++ b/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)) diff --git a/doc/stdlib/base.rst b/doc/stdlib/base.rst index fe9ffce0c78b52..81516b7a998dff 100644 --- a/doc/stdlib/base.rst +++ b/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. diff --git a/doc/stdlib/linalg.rst b/doc/stdlib/linalg.rst index 39c9dd92164142..bb5b0f1af8890b 100644 --- a/doc/stdlib/linalg.rst +++ b/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 diff --git a/doc/stdlib/pkg.rst b/doc/stdlib/pkg.rst index 492b8c27fb519f..c1f6ed83cc6332 100644 --- a/doc/stdlib/pkg.rst +++ b/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) diff --git a/test/linalg/cholesky.jl b/test/linalg/cholesky.jl index 76725cfc5e133e..8596b0f3c54e2c 100644 --- a/test/linalg/cholesky.jl +++ b/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