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
IntervalMatrix	codac2::infinite_sum_enclosure (const IntervalMatrix &A, double &mrad)
	Compute an upper bound of A+A^2+A^3+..., with A a matrix of intervals as an "error term" (use only bounds on coefficients)
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 approximate inverse of a matrix.
template<typename OtherDerived>
IntervalMatrix	codac2::inverse_enclosure (const Eigen::MatrixBase< OtherDerived > &A)
	Enclosure of the inverse of a (non-singular) matrix expression.
IntervalMatrix	codac2::inverse_enclosure (const IntervalMatrix &A)
	Enclosure of the inverse of a matrix of intervals.

Detailed Description

Date

2024

Author

Damien Massé

Copyright

Copyright 2024 Codac Team

License: GNU Lesser General Public License (LGPL)

Function Documentation

infinite_sum_enclosure()

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

Compute an upper bound of A+A^2+A^3+..., with A a matrix of intervals as an "error term" (use 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 = oo, then the inversion calculation (e.g., performed by Eigen) has somehow failed and some coefficients of the output interval matrix are [-oo,+oo]).

Precondition

A is a square matrix

Parameters

A	a matrix of intervals (supposed around 0)
mrad	the maximum radius of the result added (output argument)

Returns

the enclosure. May include (-oo,oo)

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 approximate inverse of a matrix.

Precondition

A and B are square matrices

Template Parameters

O	if LEFT_INV, use the inverse of BA (otherwise use the inverse of AB, left inverse is normally better)

Parameters

A	a matrix expression
B	a (almost punctual) approximation of its inverse,

Returns

the enclosure

  {
    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() [1/2]

template<typename OtherDerived>

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

inline

Enclosure of the inverse of a (non-singular) matrix expression.

Precondition

A is a square matrix

Parameters

A	a matrix expression

Returns

the enclosure. Can have (-oo,oo) coefficients if A is singular or almost singular

  {
    assert_release(A.is_squared());
    Index N=A.rows();
    return inverse_correction<LEFT_INV>(A, 
      A.fullPivLu().solve(Matrix::Identity(N,N)));
  }

inverse_enclosure() [2/2]

IntervalMatrix codac2::inverse_enclosure

(

const IntervalMatrix &

A

)

Enclosure of the inverse of a matrix of intervals.

Precondition

A is a square matrix

Parameters

A	a matrix of intervals

Returns

the enclosure. Can have (-oo,oo) coefficients if the inversion "failed"