Reliable inversions of matrices

Main author: Damien Massé

template<typename OtherDerived> inline IntervalMatrix codac2::inverse_enclosure(const Eigen::MatrixBase<OtherDerived> &A)

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

Parameters:: A – A matrix expression, possibly interval.
Pre:: \(\mathbf{A}\) is a square matrix.
Returns:: The enclosure of the inverse. Can have \((-\infty,\infty)\) coefficients if \(\mathbf{A}\) is singular or almost singular, if the inversion “failed”.

A = Matrix([
  [ 1, 2, 0 ],
  [ 3, 4, 1 ],
  [ 0, 1, 0 ],
])

B = inverse_enclosure(A)
# B == [[     <1, 1> ,    <-0, 0> ,        <-2, -2> ]
#       [    <-0, 0> ,    <-0, 0> , [0.9999, 1.001] ]
#       [   <-3, -3> ,     <1, 1> ,  [1.999, 2.001] ]]

i = (A*B).contains(Matrix.eye(3,3))
# i == True

Matrix A({
  { 1, 2, 0 },
  { 3, 4, 1 },
  { 0, 1, 0 },
});

IntervalMatrix B = inverse_enclosure(A);
// B == [[     <1, 1> ,    <-0, 0> ,        <-2, -2> ]
//       [    <-0, 0> ,    <-0, 0> , [0.9999, 1.001] ]
//       [   <-3, -3> ,     <1, 1> ,  [1.999, 2.001] ]]

bool i = (A.template cast<Interval>()*B).contains(Matrix::eye(3,3));
// i == true

A = Matrix({ ...
  { 1, 2, 0 }, ...
  { 3, 4, 1 }, ...
  { 0, 1, 0 }, ...
});

B = inverse_enclosure(A);
% B == [[     <1, 1> ,    <-0, 0> ,        <-2, -2> ]
%       [    <-0, 0> ,    <-0, 0> , [0.9999, 1.001] ]
%       [   <-3, -3> ,     <1, 1> ,  [1.999, 2.001] ]]

C = A*B;
i = C.contains(Matrix(3,3).eye(3,3));
% i == 1 (true)

template<LeftOrRightInv O = LEFT_INV, typename OtherDerived, typename OtherDerived_> inline 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 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.

Pre:

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

Returns:

The enclosure of the inverse.

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]\).

Parameters:

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

Pre:

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

Returns:

The sum enclosure. May be unbounded.