From 5a98276386594e13598ba501ade9adb5d8998eec Mon Sep 17 00:00:00 2001 From: damien-masse Date: Wed, 30 Sep 2026 10:40:59 +0200 Subject: [PATCH] Adding determinant operator (fwd and centered version) for dimension >=3. Adding cofactor matrix computation. Adding bwd x.M operator (x=scalar, M=matrix). --- .../functions/analytic/analytic_operators.rst | 5 +- .../src/core/matrices/codac2_py_inversion.cpp | 8 + src/core/matrices/codac2_inversion.h | 147 ++++++++++++++++++ src/core/operators/codac2_arith_mul.h | 8 +- src/core/operators/codac2_det.h | 62 +++++--- .../core/matrices/codac2_tests_inversion.cpp | 7 + tests/core/matrices/codac2_tests_inversion.py | 16 +- 7 files changed, 228 insertions(+), 25 deletions(-) diff --git a/doc/manual/manual/functions/analytic/analytic_operators.rst b/doc/manual/manual/functions/analytic/analytic_operators.rst index f1c8a3a19..2539a344f 100644 --- a/doc/manual/manual/functions/analytic/analytic_operators.rst +++ b/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. diff --git a/python/src/core/matrices/codac2_py_inversion.cpp b/python/src/core/matrices/codac2_py_inversion.cpp index 8c513e7e2..20b716e8a 100644 --- a/python/src/core/matrices/codac2_py_inversion.cpp +++ b/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) + ; } diff --git a/src/core/matrices/codac2_inversion.h b/src/core/matrices/codac2_inversion.h index 9f014c7b9..94f3df0a9 100644 --- a/src/core/matrices/codac2_inversion.h +++ b/src/core/matrices/codac2_inversion.h @@ -10,8 +10,10 @@ #pragma once #include +#include #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 + inline std::pair + cofactor_matrix_enclosure(const Eigen::MatrixBase& A) { + assert_release(A.is_squared()); + + Index N = A.cols(); + auto A_ = A.template cast(); + + Eigen::FullPivLU 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=0;r--) { + Interval s(mLU(r,c)); + s/=mLU(c,c); + for (int k=r+1;k() = ImLU; + IntervalMatrix InvU = + IntervalMatrix::Identity(N,N); + InvU.triangularView() = 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 =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;l2) { + /* build the cofactor matrix */ + Interval eps=0.0; + for (Index c=lastcol;ceps.ub()) eps=nrm; + } + Interval v1 = pow(eps,(int)a)*Interval(-1.0,1.0); + Interval d1 = v1*eps; + for (Index l=0;l2"); +// 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); @@ -110,8 +111,10 @@ namespace codac2 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().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 }; } diff --git a/tests/core/matrices/codac2_tests_inversion.cpp b/tests/core/matrices/codac2_tests_inversion.cpp index d5b553d49..1592bdef4 100644 --- a/tests/core/matrices/codac2_tests_inversion.cpp +++ b/tests/core/matrices/codac2_tests_inversion.cpp @@ -51,4 +51,11 @@ TEST_CASE("Matrix") y = inverse_enclosure(w); CHECK((w.template cast()*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))); + + } diff --git a/tests/core/matrices/codac2_tests_inversion.py b/tests/core/matrices/codac2_tests_inversion.py index e4bf2b68c..2661d1658 100644 --- a/tests/core/matrices/codac2_tests_inversion.py +++ b/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() \ No newline at end of file + unittest.main()