diff --git a/lib/src/Base/Type/MatrixImplementation.cxx b/lib/src/Base/Type/MatrixImplementation.cxx index b4b3095633..6e984fa857 100644 --- a/lib/src/Base/Type/MatrixImplementation.cxx +++ b/lib/src/Base/Type/MatrixImplementation.cxx @@ -895,6 +895,30 @@ MatrixImplementation MatrixImplementation::solveLinearSystemSquare(const MatrixI return A.solveLinearSystemSquareInPlace(b); } +/* Square matrix inverse */ +MatrixImplementation MatrixImplementation::inverseSquare() const +{ + if (nbColumns_ != nbRows_ ) throw InvalidDimensionException(HERE) << "The matrix has " << nbRows_ << " and " << nbColumns_ << ", expected a square matrix."; + MatrixImplementation identity(nbRows_, nbColumns_); + for(UnsignedInteger i = 0; i < nbRows_; ++i) + identity(i, i) = 1.0; + MatrixImplementation A(*this); + const MatrixImplementation inverseMatrix(A.solveLinearSystemSquareInPlace(identity)); + return inverseMatrix; +} + +/* Symmetric matrix inverse */ +MatrixImplementation MatrixImplementation::inverseSym() const +{ + if (nbColumns_ != nbRows_ ) throw InvalidDimensionException(HERE) << "The matrix has " << nbRows_ << " and " << nbColumns_ << ", expected a square matrix."; + MatrixImplementation identity(nbRows_, nbColumns_); + for(UnsignedInteger i = 0; i < nbRows_; ++i) + identity(i, i) = 1.0; + MatrixImplementation A(*this); + const MatrixImplementation inverseMatrix(A.solveLinearSystemSymInPlace(identity)); + return inverseMatrix; +} + /* Resolution of a linear system : square matrix */ Point MatrixImplementation::solveLinearSystemSquareInPlace(const Point & b) { diff --git a/lib/src/Base/Type/SquareMatrix.cxx b/lib/src/Base/Type/SquareMatrix.cxx index 11ae4a69ab..09177dad46 100644 --- a/lib/src/Base/Type/SquareMatrix.cxx +++ b/lib/src/Base/Type/SquareMatrix.cxx @@ -264,4 +264,11 @@ Bool SquareMatrix::isDiagonal() const return true; } + +/** Compute inverse */ +SquareMatrix SquareMatrix::inverse() const +{ + return getImplementation()->inverseSquare(); +} + END_NAMESPACE_OPENTURNS diff --git a/lib/src/Base/Type/SymmetricMatrix.cxx b/lib/src/Base/Type/SymmetricMatrix.cxx index b9556f96d5..3e5f1b5545 100644 --- a/lib/src/Base/Type/SymmetricMatrix.cxx +++ b/lib/src/Base/Type/SymmetricMatrix.cxx @@ -348,5 +348,10 @@ Scalar SymmetricMatrix::computeSumElements() const return getImplementation()->computeSumElements(); } +/** Compute inverse */ +SymmetricMatrix SymmetricMatrix::inverse() const +{ + return getImplementation()->inverseSym(); +} END_NAMESPACE_OPENTURNS diff --git a/lib/src/Base/Type/openturns/MatrixImplementation.hxx b/lib/src/Base/Type/openturns/MatrixImplementation.hxx index 5a83fee8c7..9f437b3658 100644 --- a/lib/src/Base/Type/openturns/MatrixImplementation.hxx +++ b/lib/src/Base/Type/openturns/MatrixImplementation.hxx @@ -202,6 +202,12 @@ public: MatrixImplementation solveLinearSystemSquareInPlace(const MatrixImplementation & b); MatrixImplementation solveLinearSystemSquare(const MatrixImplementation & b) const; + /** Square inverse */ + MatrixImplementation inverseSquare() const; + + /** Symmetric inverse */ + MatrixImplementation inverseSym() const; + /** Resolution of a linear system in case of a triangular matrix */ Point solveLinearSystemTri(const Point & b, const Bool lower = true, diff --git a/lib/src/Base/Type/openturns/SquareMatrix.hxx b/lib/src/Base/Type/openturns/SquareMatrix.hxx index 8c3a4dd9dc..eac2a50eac 100644 --- a/lib/src/Base/Type/openturns/SquareMatrix.hxx +++ b/lib/src/Base/Type/openturns/SquareMatrix.hxx @@ -1,6 +1,6 @@ // -*- C++ -*- /** - * @brief SquareMatrix implements the classical mathematical square matrix + * @brief SquareMatrix implements the square matrix * * Copyright 2005-2024 Airbus-EDF-IMACS-ONERA-Phimeca * @@ -35,7 +35,7 @@ class SquareComplexMatrix; /** * @class SquareMatrix * - * SquareMatrix implements the classical mathematical square matrix + * SquareMatrix implements the square matrix */ class OT_API SquareMatrix : @@ -150,6 +150,9 @@ public: /** Check if it is diagonal */ Bool isDiagonal() const; + /** Inverse matrix*/ + SquareMatrix inverse() const; + protected: private: diff --git a/lib/src/Base/Type/openturns/SymmetricMatrix.hxx b/lib/src/Base/Type/openturns/SymmetricMatrix.hxx index 5d420b4a32..2e16cc2d3d 100644 --- a/lib/src/Base/Type/openturns/SymmetricMatrix.hxx +++ b/lib/src/Base/Type/openturns/SymmetricMatrix.hxx @@ -1,6 +1,6 @@ // -*- C++ -*- /** - * @brief SymmetricMatrix implements the classical mathematical symmetric matrix + * @brief SymmetricMatrix implements the symmetric matrix * * Copyright 2005-2024 Airbus-EDF-IMACS-ONERA-Phimeca * @@ -31,7 +31,7 @@ class IdentityMatrix; /** * @class SymmetricMatrix * - * SymmetricMatrix implements the classical mathematical square matrix + * SymmetricMatrix implements the symmetric matrix */ class OT_API SymmetricMatrix : @@ -169,6 +169,9 @@ public: /** Sum all coefficients */ Scalar computeSumElements() const override; + /** Inverse matrix*/ + SymmetricMatrix inverse() const; + protected: diff --git a/lib/test/t_SquareMatrix_std.cxx b/lib/test/t_SquareMatrix_std.cxx index b88563d996..854e80e943 100644 --- a/lib/test/t_SquareMatrix_std.cxx +++ b/lib/test/t_SquareMatrix_std.cxx @@ -169,6 +169,29 @@ int main(int, char *[]) << "squareMatrix1 is empty = " << squareMatrix1.isEmpty() << std::endl << "squareMatrix5 is empty = " << squareMatrix5.isEmpty() << std::endl; + /* Check inverse() */ + SquareMatrix squareMatrix6(3); + squareMatrix6(0, 0) = 1.0; + squareMatrix6(0, 1) = 2.0; + squareMatrix6(0, 2) = 3.0; + squareMatrix6(1, 0) = 3.0; + squareMatrix6(1, 1) = 2.0; + squareMatrix6(1, 2) = 1.0; + squareMatrix6(2, 0) = 2.0; + squareMatrix6(2, 1) = 1.0; + squareMatrix6(2, 2) = 3.0; + SquareMatrix squareMatrix7(squareMatrix6.inverse()); + SquareMatrix inverseReference(3); + inverseReference(0, 0) = -5.0 / 12.0; + inverseReference(0, 1) = 3.0 / 12.0; + inverseReference(0, 2) = 4.0 / 12.0; + inverseReference(1, 0) = 7.0 / 12.0; + inverseReference(1, 1) = 3.0 / 12.0; + inverseReference(1, 2) = -8.0 / 12.0; + inverseReference(2, 0) = 1.0 / 12.0; + inverseReference(2, 1) = -3.0 / 12.0; + inverseReference(2, 2) = 4.0 / 12.0; + assert_almost_equal(squareMatrix7, inverseReference); } catch (TestFailed & ex) { diff --git a/lib/test/t_SymmetricMatrix_std.cxx b/lib/test/t_SymmetricMatrix_std.cxx index e2fc12bc1a..b9eb0fb3ed 100644 --- a/lib/test/t_SymmetricMatrix_std.cxx +++ b/lib/test/t_SymmetricMatrix_std.cxx @@ -220,6 +220,30 @@ int main(int, char *[]) M3(0, 2) = 7.1; M3(1, 2) = 9.0; fullprint << "SM * M3 = " << SM * M3 << std::endl; + + /* Check inverse() */ + SymmetricMatrix symMatrix6(3); + symMatrix6(0, 0) = 4.0; + symMatrix6(0, 1) = 2.0; + symMatrix6(0, 2) = 1.0; + symMatrix6(1, 0) = 2.0; + symMatrix6(1, 1) = 5.0; + symMatrix6(1, 2) = 3.0; + symMatrix6(2, 0) = 1.0; + symMatrix6(2, 1) = 3.0; + symMatrix6(2, 2) = 6.0; + SymmetricMatrix symMatrix7(symMatrix6.inverse()); + SymmetricMatrix inverseReference(3); + inverseReference(0, 0) = 21.0 / 67.0; + inverseReference(0, 1) = -9.0 / 67.0; + inverseReference(0, 2) = 1.0 / 67.0; + inverseReference(1, 0) = -9.0 / 67.0; + inverseReference(1, 1) = 23.0 / 67.0; + inverseReference(1, 2) = -10.0 / 67.0; + inverseReference(2, 0) = 1.0 / 67.0; + inverseReference(2, 1) = -10.0 / 67.0; + inverseReference(2, 2) = 16.0 / 67.0; + assert_almost_equal(symMatrix7, inverseReference); } catch (TestFailed & ex) { diff --git a/python/src/Matrix_doc.i.in b/python/src/Matrix_doc.i.in index a66f61437e..1c8e13369c 100644 --- a/python/src/Matrix_doc.i.in +++ b/python/src/Matrix_doc.i.in @@ -82,14 +82,14 @@ The economic QR factorization of a rectangular matrix :math:`\mat{M}` with .. math:: \mat{M} = \mat{Q} \mat{R} - = \mat{Q} \begin{bmatrix} \mat{R_1} \\ \mat{0} \end{bmatrix} - = \begin{bmatrix} \mat{Q_1}, \mat{Q_2} \end{bmatrix} - \begin{bmatrix} \mat{R_1} \\ \mat{0} \end{bmatrix} - = \mat{Q_1} \mat{R_1} - -where :math:`\mat{R_1}` is an :math:`n_c \times n_c` upper triangular matrix, -:math:`\mat{Q_1}` is :math:`n_r \times n_c`, :math:`\mat{Q_2}` is -:math:`n_r \times (n_r - n_c)`, and :math:`\mat{Q_1}` and :math:`\mat{Q_2}` + = \mat{Q} \begin{bmatrix} \mat{R}_1 \\ \mat{0} \end{bmatrix} + = \begin{bmatrix} \mat{Q}_1, \mat{Q}_2 \end{bmatrix} + \begin{bmatrix} \mat{R}_1 \\ \mat{0} \end{bmatrix} + = \mat{Q}_1 \mat{R}_1 + +where :math:`\mat{R}_1` is an :math:`n_c \times n_c` upper triangular matrix, +:math:`\mat{Q}_1` is :math:`n_r \times n_c`, :math:`\mat{Q}_2` is +:math:`n_r \times (n_r - n_c)`, and :math:`\mat{Q}_1` and :math:`\mat{Q}_2` both have orthogonal columns. Parameters @@ -118,7 +118,7 @@ least-square sense because it implies solving a (simple) triangular system: .. math:: \vect{\hat{x}} = \arg\min\limits_{\vect{x} \in \Rset^{n_r}} \|\mat{M} \vect{x} - \vect{b}\| - = \mat{R_1}^{-1} (\Tr{\mat{Q_1}} \vect{b}) + = \mat{R}_1^{-1} (\Tr{\mat{Q}_1} \vect{b}) This uses LAPACK's `DGEQRF `_ and `DORGQR `_. @@ -423,8 +423,8 @@ MMT : :class:`~openturns.Matrix` Notes ----- -When transposed is set to `True`, the method computes :math:`cM^t \times \cM`. -Otherwise it computes :math:`\cM \ times \cM^t` +When `transposed` is `True`, compute :math:`\Tr{M} M`. +Otherwise, compute :math:`M \Tr{M}`. Examples -------- @@ -475,12 +475,15 @@ C : :class:`~openturns.Matrix` Notes ----- -The matrix :math:`\cC` resulting from the Hadamard product of the matrices -:math:`\cA` and :`\cB` is as follows: +The matrix :math:`\mat{C} \in \Rset^{m \times n}` resulting from the Hadamard product ( +also known as the elementwise product) of the matrices +:math:`\mat{A} \in \Rset^{m \times n}` and :math:`\mat{B} \in \Rset^{m \times n}` is: .. math:: - \cC_{i,j} = \cA_{i,j} * \cB_{i,j} + c_{i,j} = a_{i,j} b_{i,j} + +for any :math:`i = 1, \cdots, m` and :math:`j = 1, \cdots, n`. Examples -------- @@ -508,11 +511,11 @@ sum : a float Notes ----- -We compute here the sum of elements of the matrix :math:`\cM`, that defines as: +Compute the sum of elements of the matrix :math:`\mat{A} \in \Rset^{m \times n}`: .. math:: - s = \sum_{i=1}^{nRows}\sum_{j=1}^{nColumns} \cM_{i,j} + s = \sum_{i=1}^{m} \sum_{j=1}^{n} a_{i,j} Examples -------- diff --git a/python/src/SquareMatrix_doc.i.in b/python/src/SquareMatrix_doc.i.in index 7152a704b9..c0333ab041 100644 --- a/python/src/SquareMatrix_doc.i.in +++ b/python/src/SquareMatrix_doc.i.in @@ -32,7 +32,7 @@ Get or set terms [[ 1 ] [ 1 ]] -Create an openturns matrix from a **square** numpy 2d-array (or matrix, or +Create a matrix from a **square** Numpy 2d-array (or matrix, or 2d-list)... >>> import numpy as np @@ -308,3 +308,22 @@ Examples >>> M = ot.SquareMatrix([[1.0, 2.0], [3.0, 4.0]]) >>> M.computeTrace() 5.0" + +// --------------------------------------------------------------------- + +%feature("docstring") OT::SquareMatrix::inverse +"Compute the inverse of the matrix. + +Returns +------- +inverseMatrix : :class:`~openturns.SquareMatrix` + The inverse of the matrix. + +Examples +-------- +>>> import openturns as ot +>>> M = ot.SquareMatrix([[1.0, 2.0, 3.0], [3.0, 2.0, 1.0], [2.0, 1.0, 3.0]]) +>>> print(12.0 * M.inverse()) +[[ -5 3 4 ] + [ 7 3 -8 ] + [ 1 -3 4 ]]" diff --git a/python/src/SymmetricMatrix_doc.i.in b/python/src/SymmetricMatrix_doc.i.in index 29004a1497..ce9941818e 100644 --- a/python/src/SymmetricMatrix_doc.i.in +++ b/python/src/SymmetricMatrix_doc.i.in @@ -36,7 +36,7 @@ Get or set terms [[ 1 ] [ 2 ]] -Create an openturns matrix from a **symmetric** numpy 2d-array (or matrix, or +Create a matrix from a **symmetric** Numpy 2d-array (or matrix, or 2d-list)... >>> import numpy as np @@ -158,3 +158,22 @@ Examples %feature("docstring") OT::SymmetricMatrix::checkSymmetry "Check if the internal representation is really symmetric." + +// --------------------------------------------------------------------- + +%feature("docstring") OT::SymmetricMatrix::inverse +"Compute the inverse of the matrix. + +Returns +------- +inverseMatrix : :class:`~openturns.SymmetricMatrix` + The inverse of the matrix. + +Examples +-------- +>>> import openturns as ot +>>> M = ot.SymmetricMatrix([[4.0, 2.0, 1.0], [2.0, 5.0, 3.0], [1.0, 3.0, 6.0]]) +>>> print(67.0 * M.inverse()) +[[ 21 -9 1 ] + [ -9 23 -10 ] + [ 1 -10 16 ]]" diff --git a/python/test/t_SquareMatrix_std.py b/python/test/t_SquareMatrix_std.py index ee599decdf..a3306258cc 100755 --- a/python/test/t_SquareMatrix_std.py +++ b/python/test/t_SquareMatrix_std.py @@ -1,6 +1,7 @@ #! /usr/bin/env python import openturns as ot +from openturns.testing import assert_almost_equal ot.TESTPREAMBLE() @@ -126,3 +127,12 @@ print("squareMatrix0 is empty = ", squareMatrix0.isEmpty()) print("squareMatrix1 is empty = ", squareMatrix1.isEmpty()) print("squareMatrix5 is empty = ", squareMatrix5.isEmpty()) + +# Check inverse() +squareMatrix6 = ot.SquareMatrix([[1.0, 2.0, 3.0], [3.0, 2.0, 1.0], [2.0, 1.0, 3.0]]) +squareMatrix7 = squareMatrix6.inverse() +inverseReference = ot.SquareMatrix( + [[-5.0, 3.0, 4.0], [7.0, 3.0, -8.0], [1.0, -3.0, 4.0]] +) +inverseReference /= 12.0 +assert_almost_equal(squareMatrix7, inverseReference) diff --git a/python/test/t_SymmetricMatrix_std.py b/python/test/t_SymmetricMatrix_std.py index 2956ce56bd..638a3b7706 100755 --- a/python/test/t_SymmetricMatrix_std.py +++ b/python/test/t_SymmetricMatrix_std.py @@ -1,6 +1,7 @@ #! /usr/bin/env python import openturns as ot +from openturns.testing import assert_almost_equal ot.TESTPREAMBLE() @@ -127,3 +128,10 @@ print("symmetricMatrix0 is empty = ", symmetricMatrix0.isEmpty()) print("symmetricMatrix1 is empty = ", symmetricMatrix1.isEmpty()) print("symmetricMatrix5 is empty = ", symmetricMatrix5.isEmpty()) + +# Check inverse() +symmetricMatrix6 = ot.SymmetricMatrix([[4.0, 2.0, 1.0], [2.0, 5.0, 3.0], [1.0, 3.0, 6.0]]) +symmetricMatrix7 = symmetricMatrix6.inverse() +inverseReference = ot.SymmetricMatrix([[21.0, -9.0, 1.0], [-9.0, 23.0, -10.0], [1.0, -10.0, 16.0]]) +inverseReference /= 67.0 +assert_almost_equal(symmetricMatrix7, inverseReference)