/**
* 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<number> &v);
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);
+ }
+ }
+
}
}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/lac/lapack_full_matrix.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <iostream>
+
+
+template <typename NumberType>
+void
+test()
+{
+ const unsigned int size = 4;
+ LAPACKFullMatrix<NumberType> A(size);
+ A.set_property(LAPACKSupport::symmetric);
+ Vector<NumberType> 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<double>();
+}
--- /dev/null
+
+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
projects / dealii.git / commitdiff