]> https://gitweb.dealii.org/ - dealii.git/commitdiff
implement downdating of Cholesky factorization
authorDenis Davydov <davydden@gmail.com>
Mon, 18 Dec 2017 15:52:15 +0000 (16:52 +0100)
committerDenis Davydov <davydden@gmail.com>
Tue, 19 Dec 2017 12:27:14 +0000 (13:27 +0100)
include/deal.II/lac/lapack_full_matrix.h
source/lac/lapack_full_matrix.cc
tests/lapack/full_matrix_29.cc [new file with mode: 0644]
tests/lapack/full_matrix_29.output [new file with mode: 0644]

index f63a7dae6881bc734bdb0af1812a4fbba37ce5a9..0479b3ba0c11ad65a5fe67cb33a9f252b88b7440 100644 (file)
@@ -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<number> &v);
index c2007fd7cb08ba4c3b8151b36ded79b27f3eba3f..77d18620bad6cf2797141be8a80c229d4277d260 100644 (file)
@@ -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 (file)
index 0000000..f7569b6
--- /dev/null
@@ -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 <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>();
+}
diff --git a/tests/lapack/full_matrix_29.output b/tests/lapack/full_matrix_29.output
new file mode 100644 (file)
index 0000000..77c4a5d
--- /dev/null
@@ -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  

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.