codac2_inversion.h File Reference

#include <ostream>
#include "codac2_Matrix.h"
#include "codac2_IntervalMatrix.h"

Include dependency graph for codac2_inversion.h:

Go to the source code of this file.

Functions
template<LeftOrRightInv O = LEFT_INV, typename OtherDerived, typename OtherDerived_>
IntervalMatrix	codac2::inverse_correction (const Eigen::MatrixBase< OtherDerived > &A, const Eigen::MatrixBase< OtherDerived_ > &B)
	Correct the approximation of the inverse \(\mathbf{B}\approx\mathbf{A}^{-1}\) of a square matrix \(\mathbf{A}\) by providing a reliable enclosure \([\mathbf{A}^{-1}]\).
template<typename OtherDerived>
IntervalMatrix	codac2::inverse_enclosure (const Eigen::MatrixBase< OtherDerived > &A)
	Enclosure of the inverse of a (non-singular) matrix expression, possibly an interval matrix.
IntervalMatrix	codac2::infinite_sum_enclosure (const IntervalMatrix &A, double &mrad)
	Compute an upper bound of \(\left([\mathbf{A}]+[\mathbf{A}]^2+[\mathbf{A}]^3+\dots\right)\), with \([\mathbf{A}]\) a matrix of intervals as an "error term" (uses only bounds on coefficients).

Detailed Description

Date

2024

Author

Damien Massé

Copyright

Copyright 2024 Codac Team

License: GNU Lesser General Public License (LGPL)

Function Documentation

inverse_correction()

template<LeftOrRightInv O = LEFT_INV, typename OtherDerived, typename OtherDerived_>

IntervalMatrix codac2::inverse_correction ( const Eigen::MatrixBase< OtherDerived > & A, const Eigen::MatrixBase< OtherDerived_ > & B )

inline

Correct the approximation of the inverse \(\mathbf{B}\approx\mathbf{A}^{-1}\) of a square matrix \(\mathbf{A}\) by providing a reliable enclosure \([\mathbf{A}^{-1}]\).

Precondition

\(\mathbf{A}\) and \(\mathbf{B}\) are square matrices, possibly interval matrices.

Template Parameters

O	If LEFT_INV, use the inverse of \(\mathbf{BA}\) (otherwise use the inverse of \(\mathbf{AB}\), left inverse is normally better). In Python/Matlab, this template parameter is provided as a last boolean argument, and is left_inv = True by default.

Parameters

A	A matrix expression, possibly interval.
B	An (almost punctual) approximation of its inverse.

Returns

The enclosure of the inverse.

  {
    assert_release(A.is_squared());
    assert_release(B.is_squared());
 
    auto A_ = A.template cast<Interval>();
    auto B_ = B.template cast<Interval>();
 
    Index N = A_.rows();
    assert_release(N==B_.rows());
 
    auto Id = IntervalMatrix::Identity(N,N);
    auto erMat = [&]() { if constexpr(O == LEFT_INV) return -B_*A_+Id; else return -A_*B_+Id; }();
    
    double mrad=0.0;
    IntervalMatrix E = infinite_sum_enclosure(erMat,mrad);
    IntervalMatrix Ep = Id+erMat*(Id+E); 
        /* one could use several iterations here, either
           using mrad, or directly */
 
    auto res = (O == LEFT_INV) ? IntervalMatrix(Ep*B_) : IntervalMatrix(B_*Ep);
 
    // We want the presence of non-invertible matrices to
    // result in coefficients of the form [-oo,+oo].
    if (mrad==oo) {
       for (Index c=0;c<N;c++) {
         for (Index r=0;r<N;r++) {
           if (Ep(r,c).is_unbounded()) {
              for (Index k=0;k<N;k++) {
                   if constexpr(O == LEFT_INV)
                       res(r,k) = Interval();
                   else 
                       res(k,c) = Interval();
              }
           }
         }
       }
    }
    return  res;
  }

inverse_enclosure()

template<typename OtherDerived>

IntervalMatrix codac2::inverse_enclosure ( const Eigen::MatrixBase< OtherDerived > & A )

inline

Enclosure of the inverse of a (non-singular) matrix expression, possibly an interval matrix.

Precondition

\(\mathbf{A}\) is a square matrix.

Parameters

A	A matrix expression, possibly interval.

Returns

The enclosure of the inverse. Can have \((-\infty,\infty)\) coefficients if \(\mathbf{A}\) is singular or almost singular, if the inversion "failed".

  {
    assert_release(A.is_squared());
    Index N = A.rows();
 
    if constexpr(std::is_same_v<typename OtherDerived::Scalar,Interval>)
      return inverse_correction<LEFT_INV>(A, 
        (A.mid()).fullPivLu().solve(Matrix::Identity(N,N)));
 
    else
      return inverse_correction<LEFT_INV>(A, 
        A.fullPivLu().solve(Matrix::Identity(N,N)));
  }

infinite_sum_enclosure()

IntervalMatrix codac2::infinite_sum_enclosure	(	const IntervalMatrix &	A,
		double &	mrad )

Compute an upper bound of \(\left([\mathbf{A}]+[\mathbf{A}]^2+[\mathbf{A}]^3+\dots\right)\), with \([\mathbf{A}]\) a matrix of intervals as an "error term" (uses only bounds on coefficients).

The function also returns mrad, which gives an idea of the magnification of the matrix during calculation. In particular, if mrad = \(\infty\), then the inversion calculation (e.g., performed by Eigen) has somehow failed and some coefficients of the output interval matrix are \([-\infty,\infty]\).

Precondition

\([\mathbf{A}]\) is a square matrix.

Parameters

A	A matrix of intervals (supposed around \(\mathbf{0}\)).
mrad	The maximum radius of the result added (output argument).

Returns

The sum enclosure. May be unbounded.