From b965a781b940765c3b62a4fd00adfeabe5bc18f2 Mon Sep 17 00:00:00 2001 From: Denis Davydov Date: Mon, 18 Dec 2017 16:52:15 +0100 Subject: [PATCH] implement downdating of Cholesky factorization --- include/deal.II/lac/lapack_full_matrix.h | 5 + source/lac/lapack_full_matrix.cc | 50 ++++++++- tests/lapack/full_matrix_29.cc | 129 +++++++++++++++++++++++ tests/lapack/full_matrix_29.output | 16 +++ 4 files changed, 196 insertions(+), 4 deletions(-) create mode 100644 tests/lapack/full_matrix_29.cc create mode 100644 tests/lapack/full_matrix_29.output diff --git a/include/deal.II/lac/lapack_full_matrix.h b/include/deal.II/lac/lapack_full_matrix.h index f63a7dae68..0479b3ba0c 100644 --- a/include/deal.II/lac/lapack_full_matrix.h +++ b/include/deal.II/lac/lapack_full_matrix.h @@ -144,6 +144,11 @@ public: /** * Perform a rank-1 update of a symmetric matrix * $ A \leftarrow A + a \, \rm v \rm v^T $. + * + * This function also works for Cholesky factorization. Updating ($a>0$) is + * performed via Givens rotations, whereas downdating ($a<0$) via hyperbolic + * rotations. Note that the later case might lead to a negative definite + * matrix in which case the error will be thrown. */ void add(const number a, const Vector &v); diff --git a/source/lac/lapack_full_matrix.cc b/source/lac/lapack_full_matrix.cc index c2007fd7cb..77d18620ba 100644 --- a/source/lac/lapack_full_matrix.cc +++ b/source/lac/lapack_full_matrix.cc @@ -248,15 +248,57 @@ namespace const number t = A(i,k); A(i,k) = csr[0] * A(i,k) + csr[1] * z(i); z(i) = -csr[1] * t + csr[0] * z(i); - - //A(i,k) = ( t + csr[1] * z(i)) / csr[0]; - //z(i) = - csr[1] * t + csr[0] * z(i) ; } } } else { - AssertThrow(false, ExcNotImplemented()); + // downdating is not always possible as we may end up with + // negative definite matrix. If it's possible, then it boils + // down to application of hyperbolic rotations. + // Observe that + // + // | L^T |T | L^T | + // A <-- | | S | | = L L^T - z z^T + // | z^T | | z^T | + // + // |In 0 | + // S := | | + // |0 -1 | + // + // We are looking for H which is S-orthogonal (HSH^T=S) and + // can restore upper-triangular factor of the factorization of A above. + // We will use Hyperbolic rotations to do the job + // + // | c -s | | x1 | | r | + // | | = | | = | |, c^2 - s^2 = 1 + // |-s c | | x2 | | 0 | + // + // which have real solution only if x2 <= x1. + // See also Linpack's http://www.netlib.org/linpack/dchdd.f and + // https://infoscience.epfl.ch/record/161468/files/cholupdate.pdf and + // "Analysis of a recursive Least Squares Hyperbolic Rotation Algorithm for Signal Processing", Alexander, Pan, Plemmons, 1988. + z *= std::sqrt(-a); + for (unsigned int k = 0; k < N; ++k) + { + // we have Cholesky factor of SPD matrix + Assert (A(k,k) != 0, ExcInternalError()); + const number tau = z(k)/A(k,k); + AssertThrow (std::abs(tau) < 1., + ExcMessage("Cholesky downdating led to negative definite matrix")); + const number u = sqrt((1.-tau)*(1.+tau)); // <-- more stable than std::sqrt(1.-tau*tau) + const number c = 1./u; + const number s = c * tau; + const number r = u * std::abs(A(k,k)); + A(k,k) = r; + for (unsigned int i = k+1; i < N; ++i) + { + const number t = A(i,k); + A(i,k) = c * A(i,k) - s * z(i); + z(i) = -s * t + c * z(i); + } + } + } } diff --git a/tests/lapack/full_matrix_29.cc b/tests/lapack/full_matrix_29.cc new file mode 100644 index 0000000000..f7569b6947 --- /dev/null +++ b/tests/lapack/full_matrix_29.cc @@ -0,0 +1,129 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// test LAPACKFullMatrix::add() for rank1 downdate of a Cholesky factorization + +/* MWE in Octave: +A = pascal(4) +A = + + 1 1 1 1 + 1 2 3 4 + 1 3 6 10 + 1 4 10 20 + +x = [3 2 1 5]'; + +R = chol(A) +R = + + 1 1 1 1 + 0 1 2 3 + 0 0 1 3 + 0 0 0 1 + +R1 = cholupdate(R,x) + +R1 = + + 3.16228 2.21359 1.26491 5.05964 + 0.00000 1.04881 2.09762 2.66970 + 0.00000 0.00000 1.00000 3.00000 + 0.00000 0.00000 0.00000 1.80907 + +x = x / 2 + +R2 = cholupdate(R1,x,'-') + +R2 = + + 2.78388 1.97566 1.16743 4.40033 + 0.00000 1.04727 2.09454 2.67978 + 0.00000 0.00000 1.00000 3.00000 + 0.00000 0.00000 0.00000 1.79050 + +*/ + +#include "../tests.h" +#include "create_matrix.h" +#include +#include +#include + +#include + + +template +void +test() +{ + const unsigned int size = 4; + LAPACKFullMatrix A(size); + A.set_property(LAPACKSupport::symmetric); + Vector v(size); + + A(0,0) = 1; + A(0,1) = 1; + A(0,2) = 1; + A(0,3) = 1; + A(1,0) = 1; + A(1,1) = 2; + A(1,2) = 3; + A(1,3) = 4; + A(2,0) = 1; + A(2,1) = 3; + A(2,2) = 6; + A(2,3) = 10; + A(3,0) = 1; + A(3,1) = 4; + A(3,2) =10; + A(3,3) = 20; + + v(0) = 3; + v(1) = 2; + v(2) = 1; + v(3) = 5; + + const NumberType a = 1.; + + A.compute_cholesky_factorization(); + deallog << "Cholesky:" << std::endl; + A.print_formatted(deallog.get_file_stream(),5); + + A.add(a,v); + + deallog << "Cholesky updated:" << std::endl; + A.print_formatted(deallog.get_file_stream(),5); + + v /= 2.; + + A.add(-1.,v); + + deallog << "Cholesky downdated:" << std::endl; + A.print_formatted(deallog.get_file_stream(),5); + +} + + +int main() +{ + const std::string logname = "output"; + std::ofstream logfile(logname.c_str()); + logfile.precision(5); + deallog.attach(logfile); + + test(); +} diff --git a/tests/lapack/full_matrix_29.output b/tests/lapack/full_matrix_29.output new file mode 100644 index 0000000000..77c4a5d9d8 --- /dev/null +++ b/tests/lapack/full_matrix_29.output @@ -0,0 +1,16 @@ + +DEAL::Cholesky: +1.00000e+00 +1.00000e+00 1.00000e+00 +1.00000e+00 2.00000e+00 1.00000e+00 +1.00000e+00 3.00000e+00 3.00000e+00 1.00000e+00 +DEAL::Cholesky updated: +3.16228e+00 +2.21359e+00 1.04881e+00 +1.26491e+00 2.09762e+00 1.00000e+00 +5.05964e+00 2.66970e+00 3.00000e+00 1.80907e+00 +DEAL::Cholesky downdated: +2.78388e+00 +1.97566e+00 1.04727e+00 +1.16743e+00 2.09454e+00 1.00000e+00 +4.40033e+00 2.67978e+00 3.00000e+00 1.79050e+00 -- 2.39.5