Determinant operator, cofactor matrix, bwd scalar.matrix operator.

doc/manual/manual/functions/analytic/analytic_operators.rst

@@ -175,7 +175,7 @@ If you notice any mathematical operators missing from the list below, feel free
   +-----------------------------------------------------+----------------------+---------------+-------------------------------------+--------+--------+-------+------------+
   | :math:`\cosh(x)`                                    | ``cosh(x)``          | |CoshOp|      | ``x``: scalar                       ||okk|   ||okk|   ||okk|  ||okk|       |
   +-----------------------------------------------------+----------------------+---------------+-------------------------------------+--------+--------+-------+------------+
-  | :math:`\det(\mathbf{A})`                            | ``det(A)``           | |DetOp|       | ``A``: matrix                       ||bok|   ||bok|   ||bok|  ||bok|       |
+  | :math:`\det(\mathbf{A})`                            | ``det(A)``           | |DetOp|       | ``A``: matrix                       ||okk|   ||okk|   ||okk|  ||okk|       |
   +-----------------------------------------------------+----------------------+---------------+-------------------------------------+--------+--------+-------+------------+
   | :math:`\exp(x)`                                     | ``exp(x)``           | |ExpOp|       | ``x``: scalar                       ||okk|   ||okk|   ||okk|  ||okk|       |
   +-----------------------------------------------------+----------------------+---------------+-------------------------------------+--------+--------+-------+------------+
@@ -223,7 +223,7 @@ If you notice any mathematical operators missing from the list below, feel free
   +-----------------------------------------------------+----------------------+               +-------------------------------------+--------+--------+-------+------------+
   | :math:`\mathbf{x}_1\cdot x_2`                       | ``x1*x2``            |               | ``x1``: vector, ``x2``: scalar      ||okk|   ||okk|   ||okk|  ||okk|       |
   +-----------------------------------------------------+----------------------+               +-------------------------------------+--------+--------+-------+------------+
-  | :math:`x_1\cdot\mathbf{X}_2`                        | ``x1*X2``            |               | ``x1``: scalar, ``X2``: matrix      ||okk|   ||okk|   ||okk|  ||nok|       |
+  | :math:`x_1\cdot\mathbf{X}_2`                        | ``x1*X2``            |               | ``x1``: scalar, ``X2``: matrix      ||okk|   ||okk|   ||okk|  ||okk|       |
   +-----------------------------------------------------+----------------------+               +-------------------------------------+--------+--------+-------+------------+
   | :math:`\mathbf{x}_1\cdot\mathbf{x}_2`               | ``x1*x2``            |               | ``x1``: row, ``x2``: vector         ||okk|   ||nok|   ||nok|  ||okk|       |
   +-----------------------------------------------------+----------------------+               +-------------------------------------+--------+--------+-------+------------+
@@ -280,7 +280,6 @@ If you notice any mathematical operators missing from the list below, feel free
   | :math:`[\mathbf{x}](t)`                             | ``x(t)``             | |TubeOp|      | ``x``: tube, ``t``: scalar          ||okk|   ||nok|   ||nok|  ||nok|       |
   +-----------------------------------------------------+----------------------+---------------+-------------------------------------+--------+--------+-------+------------+
 
-| Note: the operator :math:`\det` is only available for :math:`1\times 1` and :math:`2\times 2` matrices.
 | Note: the operator :math:`\bmod` is only available for real periods (double precision), interval periods are not yet supported.
 
 

python/src/core/matrices/codac2_py_inversion.cpp

@@ -48,5 +48,13 @@ void export_inversion(py::module& m)
       INTERVALMATRIX_INFINITE_SUM_ENCLOSURE_CONST_INTERVALMATRIX_REF_DOUBLE_REF,
       "A"_a)
 
+    .def("cofactor_matrix_enclosure", [](const Matrix& A) { return cofactor_matrix_enclosure(A); },
+      INTERVALMATRIX_INFINITE_SUM_ENCLOSURE_CONST_INTERVALMATRIX_REF_DOUBLE_REF,
+      "A"_a)
+
+    .def("cofactor_matrix_enclosure", [](const IntervalMatrix& A) { return cofactor_matrix_enclosure(A); },
+      PAIR_INTERVALMATRIXINTERVAL_COFACTOR_MATRIX_ENCLOSURE_CONST_EIGEN_MATRIXBASE_OTHERDERIVED_REF,
+      "A"_a)
+
    ;
 }

src/core/matrices/codac2_inversion.h

@@ -10,8 +10,10 @@
 #pragma once
 
 #include <ostream>
+#include <iostream>
 #include "codac2_Matrix.h"
 #include "codac2_IntervalMatrix.h"
+#include "codac2_IntervalVector.h"
 
 namespace codac2
 {
@@ -114,4 +116,149 @@ namespace codac2
    * \return The sum enclosure. May be unbounded.
    */
   IntervalMatrix infinite_sum_enclosure(const IntervalMatrix& A, double& mrad);
+
+
+  /**
+   * \brief compute the interval cofactor matrix, as well as the 
+   * determinant, of a square matrix 
+   * \param A A matrix expression, possibly interval.
+   * \return The enclosure of the cofactor matrix (first) 
+   * and determinant (second)
+   */
+  template<typename OtherDerived>
+  inline std::pair<IntervalMatrix,Interval>
+	cofactor_matrix_enclosure(const Eigen::MatrixBase<OtherDerived>& A) {
+      assert_release(A.is_squared());
+
+      Index N = A.cols();
+      auto A_ = A.template cast<Interval>();
+
+      Eigen::FullPivLU<Matrix> lu(A_.mid());
+        /* FullPiv is used to get the ``best'' pivoting strategy */
+      auto P = lu.permutationP();
+      auto Q = lu.permutationQ();
+
+      IntervalMatrix pMq = P*A_*Q;
+      /* preconditioning ? 
+         we use L and U to compute
+         L-1 pMq U-1
+         and include L-1 and U-1 in the computation of the cofactor matrix
+         */
+      Matrix mLU = lu.matrixLU();
+      for (int i=0;i<N;i++) {
+          if (std::fabs(mLU(i,i))<1e-5) { /* FIXME : value of "nonzero" ? */
+             mLU(i,i)=(mLU(i,i)<0.0 ? -1e-5 : 1e-5); 
+          }
+      }
+      IntervalMatrix ImLU 
+                = IntervalMatrix::Zero(N,N);
+      for (int c=0;c<N;c++) {
+        for (int r=c+1;r<N;r++) {
+           Interval s(mLU(r,c));
+           for (int k=c+1;k<r;k++) {
+              s += mLU(r,k)*ImLU(k,c);
+           }
+           ImLU(r,c)=-s;
+        }
+      }
+      for (int c=0;c<N;c++) {
+        int r;
+        ImLU(c,c)=Interval(1.0)/mLU(c,c);
+        r=c-1;
+        for (;r>=0;r--) {
+           Interval s(mLU(r,c));
+           s/=mLU(c,c);
+           for (int k=r+1;k<c;k++) {
+               s += mLU(r,k)*ImLU(k,c);
+           }
+           ImLU(r,c)=-s/mLU(r,r);
+        }
+      }
+      IntervalMatrix InvL = 
+                IntervalMatrix::Identity(N,N);
+      InvL.triangularView<Eigen::StrictlyLower>() = ImLU;
+      IntervalMatrix InvU = 
+                IntervalMatrix::Identity(N,N);
+      InvU.triangularView<Eigen::Upper>() = ImLU;
+      pMq = InvL * pMq * InvU;
+
+
+
+      IntervalMatrix Linv = IntervalMatrix::Identity(N,N);
+      IntervalMatrix Uinv = IntervalMatrix::Identity(N,N);
+      /* partial inversion of pMq is done, i.e. 
+        we build a such that Linv * pMq * Uinv is equal
+         ( D 0  )   ( Id     0    )
+         ( 0 Id ) * ( 0  <eps ). The comatrix of the second part is :
+        is      ( d 0 )
+                ( 0 w ) with d = (max nrm)^a   (a = size of eps block)
+                             w = (max nrm)^(a-1) */
+      Index c=0;
+      while (c<N-1) {
+         Interval pivot = pMq(c,c);
+         if (pivot.contains(0)) break; /* or if the mig is too low */
+         Linv.block(c+1,0,N-1-c,c+1) -= 
+         pMq.block(c+1,c,N-1-c,1)*Linv.block(c,0,1,c+1)/pivot;
+         Uinv.block(0,c+1,c+1,N-1-c) -= 
+             Uinv.block(0,c,c+1,1)*pMq.block(c,c+1,1,N-1-c)/pivot;
+         pMq.bottomRightCorner(N-1-c,N-1-c) -=
+            pMq.block(c+1,c,N-1-c,1)*pMq.block(c,c+1,1,N-1-c)/pivot;
+         pMq.block(c+1,c,N-1-c,1).fill(0.0);
+         pMq.block(c,c+1,1,N-1-c).fill(0.0);
+         c++;
+      }
+      Index lastcol=c;
+      Interval det=1.0;
+      /* build the comatrix for the diagonal part of pMq 
+         (note: built in place) */
+      if (lastcol==N-1) lastcol=N; /*  the inversion succeeded completely, 
+                    pMq is diagonal  (but the last coef may contains 0) */
+      Index a = N-lastcol;
+      IntervalVector tmp = IntervalVector::Zero(lastcol);
+      tmp[0]=1.0;
+      for (Index c=0;c<lastcol-1;c++) {
+         tmp[c+1] = tmp[c]*pMq(c,c);
+      }
+      det = tmp[lastcol-1]*pMq(lastcol-1,lastcol-1);
+      for (Index c=lastcol-1;c>=1;c--) {
+         std::swap(tmp[c],pMq(c,c));
+         pMq(c,c)*=tmp[0];
+         tmp[0]*=tmp[c];
+      }
+      pMq(0,0)=tmp[0];
+      /* special case: a=2 */
+      if (a==2) {
+         Interval det2 = pMq(N-1,N-1)*pMq(N-2,N-2)-pMq(N-1,N-2)*pMq(N-2,N-1);
+         for (Index l=0;l<lastcol;l++) { 
+           pMq(l,l)*=det2;
+         }
+         pMq.bottomRightCorner(2,2) *= det;
+         std::swap(pMq(N-1,N-1),pMq(N-2,N-2));
+         std::swap(pMq(N-1,N-2),pMq(N-2,N-1));
+         pMq(N-1,N-2)=-pMq(N-1,N-2);
+         pMq(N-2,N-1)=-pMq(N-2,N-1);
+         det *= det2;
+      } else if (a>2) {
+         /* build the cofactor matrix */
+          Interval eps=0.0;
+          for (Index c=lastcol;c<N;c++) {
+             Interval nrm = pMq.block(lastcol,c,a,1).norm();
+             if (nrm.ub()>eps.ub()) eps=nrm;
+          }
+          Interval v1 = pow(eps,(int)a)*Interval(-1.0,1.0);
+          Interval d1 = v1*eps;
+          for (Index l=0;l<lastcol;l++) { 
+             pMq(l,l)*=d1;
+          }
+          pMq.bottomRightCorner(a,a).fill(v1*det);
+          det *= d1;
+      }
+      pMq = InvL.transpose()*Linv.transpose()*pMq*Uinv.transpose()*InvU.transpose();
+      return {
+        P.determinant()*Q.determinant()*
+ 	((IntervalMatrix)(P.transpose())*pMq*(IntervalMatrix(Q.transpose()))),
+        P.determinant()*Q.determinant()*det
+      };
+  }
+
 }

src/core/operators/codac2_arith_mul.h

@@ -274,9 +274,13 @@ namespace codac2
     };
   }
 
-  inline void MulOp::bwd([[maybe_unused]] const IntervalMatrix& y, [[maybe_unused]] Interval& x1, [[maybe_unused]] IntervalMatrix& x2)
+  inline void MulOp::bwd(const IntervalMatrix& y, Interval& x1, IntervalMatrix& x2)
   {
-    // todo
+    assert(y.cols() == x2.cols());
+    assert(y.rows() == x2.rows());
+    for(Index j = 0 ; j < x2.cols() ; j++)
+      for(Index i = 0 ; i < x2.rows() ; i++)
+        MulOp::bwd(y(i,j), x1, x2(i,j));
   }
 
   inline IntervalVector MulOp::fwd(const IntervalMatrix& x1, const IntervalVector& x2)

src/core/operators/codac2_det.h

@@ -16,6 +16,7 @@
 #include "codac2_AnalyticExprWrapper.h"
 #include "codac2_arith_sub.h"
 #include "codac2_arith_mul.h"
+#include "codac2_inversion.h"
 
 namespace codac2
 {
@@ -102,16 +103,18 @@ namespace codac2
   inline Interval DetOp::fwd(const IntervalMatrix& x)
   {
     assert_release(x.is_squared() && "can only compute determinants for a square matrix");
-    assert_release((x.rows() == 1 || x.rows() == 2) && "determinant not yet computable for n×n matrices, n>2");
+//    assert_release((x.rows() == 1 || x.rows() == 2) && "determinant not yet computable for n×n matrices, n>2");
 
     if(x.rows() == 1) // 1×1 matrix
       return x(0,0);
 
     else if(x.rows() == 2) // 2×2 matrix
       return x(0,0)*x(1,1)-x(0,1)*x(1,0);
 
-    else
-      return Interval::empty(); // unhandled case
+    else {
+      auto ret = cofactor_matrix_enclosure(x);
+      return ret.second;
+    }
   }
 
   inline ScalarType DetOp::fwd_natural(const MatrixType& x)
@@ -138,25 +141,35 @@ namespace codac2
 
     IntervalMatrix d(1, x.da.cols());
     if (x.a.rows()==2) { /* otherwise, will fail afterwards */
-       for (Index i=0; i < d.cols() ; i++) {
-          d(0,i) = x.da(0,i)*x.a(1,1) + x.da(3,i)*x.a(0,0)
-		 - x.da(1,i)*x.a(0,1) - x.da(2,i)*x.a(1,0);
+         for (Index i=0; i < d.cols() ; i++) {
+            d(0,i) = x.da(0,i)*x.a(1,1) + x.da(3,i)*x.a(0,0)
+  		 - x.da(1,i)*x.a(0,1) - x.da(2,i)*x.a(1,0);
 		/* note:  ColMajor is assumed here */
-       }      
+         }      
+         return {
+           fwd(x.m),
+           fwd(x.a),
+           d, 
+           x.def_domain
+         };
+    } else {
+       auto ret1 = cofactor_matrix_enclosure(x.a);
+       for (Index i=0; i < d.cols() ; i++) 
+         d(0,i)  = ret1.first.reshaped<Eigen::ColMajor>().dot(x.da.col(i));
+       return {
+         fwd(x.m),
+         ret1.second,
+         d, 
+         x.def_domain
+       };
     }
-
-    return {
-      fwd(x.m),
-      fwd(x.a),
-      d, 
-      x.def_domain
-    };
+
   }
 
   inline void DetOp::bwd(const Interval& y, IntervalMatrix& x)
   {
     assert_release(x.is_squared() && "can only compute determinants for a square matrix");
-    assert_release((x.rows() == 1 || x.rows() == 2) && "determinant not yet computable for n×n matrices, n>2");
+//    assert_release((x.rows() == 1 || x.rows() == 2) && "determinant not yet computable for n×n matrices, n>2");
 
     if(x.rows() == 1) // 1×1 matrix
       x(0,0) &= y;
@@ -169,9 +182,11 @@ namespace codac2
       MulOp::bwd(z2, x(1,0), x(0,1));
     }
 
-    else
+    else /* we only check that the determinant intersects y */
     {
-      // unhandled case
+      auto ret = cofactor_matrix_enclosure(x);
+      if (ret.second.is_disjoint(y)) x.set_empty();
+      // otherwise, do nothing 
     }
   }
 
@@ -258,10 +273,19 @@ namespace codac2
     IntervalMatrix a(3,3);
     a.col(0) = x1.a; a.col(1) = x2.a; a.col(2) = x3.a;
 
+    assert(x1.da.cols() == x2.da.cols());
+    assert(x1.da.cols() == x3.da.cols());
+    IntervalMatrix d(1, x1.da.cols());
+    auto ret1 = cofactor_matrix_enclosure(a);
+    for (Index i=0; i < d.cols() ; i++) {
+       d(0,i)  = ret1.first.col(0).dot(x1.da.col(i));
+       d(0,i)  += ret1.first.col(1).dot(x2.da.col(i));
+       d(0,i)  += ret1.first.col(2).dot(x3.da.col(i));
+    }
     return {
       fwd(m),
-      fwd(a),
-      IntervalMatrix(0,0), // not supported yet for auto diff
+      ret1.second,
+      d, 
       x1.def_domain && x2.def_domain && x3.def_domain
     };
   }

tests/core/matrices/codac2_tests_inversion.cpp

@@ -51,4 +51,11 @@ TEST_CASE("Matrix")
 
   y = inverse_enclosure(w);
   CHECK((w.template cast<Interval>()*y).contains(Matrix::Identity(3,3)));
+
+  auto cm = cofactor_matrix_enclosure(w);
+  CHECK(Approx(cm.second)==Interval(-1.0));
+  CHECK(((1.0/cm.second)*
+                cm.first.transpose()*IntervalMatrix(w)).contains(Matrix::Identity(3,3)));
+
+
 }

tests/core/matrices/codac2_tests_inversion.py

@@ -52,6 +52,20 @@ def test_Inversion_3(self):
     v = inverse_enclosure(u)
     self.assertTrue((u*v).contains(Matrix.eye(3,3)))
 
+  def test_Inversion_4(self):
+
+    w = Matrix([
+      [ 1,2,0 ],
+      [ 3,4,1 ],
+      [ 0,1,0 ]
+    ])
+
+    y = inverse_enclosure(w);
+    self.assertTrue((w*y).contains(Matrix.eye(3,3)))
+
+    cm = cofactor_matrix_enclosure(w)
+    self.assertTrue(Approx(cm[1])==-1.0)
+    self.assertTrue(((1.0/cm[1])*cm[0].transpose()*w).contains(Matrix.eye(3,3)))
 
 if __name__ ==  '__main__':
-  unittest.main()
+  unittest.main()