STK++ 0.9.13
STK_InvertSym.h
Go to the documentation of this file.
1/*--------------------------------------------------------------------*/
2/* Copyright (C) 2004-2016 Serge Iovleff, Université Lille 1, Inria
3
4 This program is free software; you can redistribute it and/or modify
5 it under the terms of the GNU Lesser General Public License as
6 published by the Free Software Foundation; either version 2 of the
7 License, or (at your option) any later version.
8
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU Lesser General Public License for more details.
13
14 You should have received a copy of the GNU Lesser General Public
15 License along with this program; if not, write to the
16 Free Software Foundation, Inc.,
17 59 Temple Place,
18 Suite 330,
19 Boston, MA 02111-1307
20 USA
21
22 Contact : S..._Dot_I..._At_stkpp_Dot_org (see copyright for ...)
23*/
24
25/*
26 * Project: stkpp::Algebra
27 * created on: 9 août 2016
28 * Author: iovleff, S..._Dot_I..._At_stkpp_Dot_org (see copyright for ...)
29 **/
30
35#ifndef STK_SYMINVERT_H
36#define STK_SYMINVERT_H
37
38#ifdef STKUSELAPACK
39#include <Algebra/include/STK_lapack_SymEigen.h>
40#else
41#include <Algebra/include/STK_SymEigen.h>
42#endif
43
44namespace STK
45{
46
47namespace hidden
48{
52template<class Matrix, int Size_>
53struct InvertSymImpl
54{
55 typedef typename Matrix::Type Type;
62 static Type invertSymMatrix11( Matrix const& m, CArraySquare<Type, Size_>& inv)
63 {
64 const int iShift = m.beginRows(), jShift = m.beginCols();
65 inv.resize(TRange<Size_>(0, 1));
66 // cofactor (0,0) [0]
67 inv(0, 0) = Type(1);
68 // compute determinant
69 Type det = m(iShift+0, jShift+0);
70 if (det == Type(0)) return Type(0);
71 // compute inverse matrix
72 inv /= det;
73 return det;
74 }
82 static Type invertSymMatrix22( Matrix const& m, CArraySquare<Type, Size_>& inv)
83 {
84 const int iShift = m.beginRows(), jShift = m.beginCols();
85 inv.resize(TRange<Size_>(0, 2));
86 // cofactor (0,0) [0]
87 inv(0, 0) = m(iShift+1, jShift+1);
88 // cofactor (1,0) [1]
89 inv(1, 0) = - m(iShift+1, jShift+0);
90 // cofactor (1,1) [3]
91 inv(1, 1) = m(iShift+0, jShift+0);
92 // symmetry
93 inv(0, 1) = inv(1,0);
94 // compute determinant
95 Type det = m(iShift+0, jShift+0) * inv(0,0)
96 + m(iShift+1, jShift+0) * inv(1,0);
97 if (det == Type(0)) return Type(0);
98 // inverse matrix
99 inv /= det;
100 return det;
101 }
108 static Type invertSymMatrix33( Matrix const& m, CArraySquare<Type, Size_>& inv)
109 {
110 const int iShift = m.beginRows(), jShift = m.beginCols();
111 inv.resize(TRange<Size_>(0, 3));
112 // cofactor (0,0) [0]
113 inv(0, 0) = m(iShift+1, jShift+1) * m(iShift+2, jShift+2)
114 - m(iShift+2, jShift+1) * m(iShift+2, jShift+1);
115 // cofactor (1,0) [1]
116 inv(1, 0) = m(iShift+1, jShift+2) * m(iShift+2, jShift+0)
117 - m(iShift+1, jShift+0) * m(iShift+2, jShift+2);
118 inv(2, 0) = m(iShift+1, jShift+0) * m(iShift+2, jShift+1)
119 - m(iShift+2, jShift+0) * m(iShift+1, jShift+1);
120 inv(1, 1) = m(iShift+0, jShift+0) * m(iShift+2, jShift+2)
121 - m(iShift+0, jShift+2) * m(iShift+2, jShift+0);
122 inv(2, 1) = m(iShift+2, jShift+0) * m(iShift+1, jShift+0)
123 - m(iShift+0, jShift+0) * m(iShift+2, jShift+1);
124 inv(2, 2) = m(iShift+0, jShift+0) * m(iShift+1, jShift+1)
125 - m(iShift+1, jShift+0) * m(iShift+0, jShift+1);
126 // symmetric
127 inv(0,1) = inv(1,0); inv(0,2) = inv(2,0);
128 inv(1,2) = inv(2,1);
129 // compute determinant
130 Type det = m(iShift+0, jShift+0) * inv(0,0)
131 + m(iShift+1, jShift+0) * inv(1,0)
132 + m(iShift+2, jShift+0) * inv(2,0);
133 if (det == Type(0)) return Type(0);
134 // inverse matrix
135 inv /= det;
136 return det;
137 }
144 static Type invertSymMatrix44( Matrix const& m, CArraySquare<Type, Size_>& inv)
145 {
146 const int iShift = m.beginRows(), jShift = m.beginCols();
147 inv.resize(TRange<Size_>(0, 4));
148 // cofactor (0,0) [0]
149 inv(0,0) = m(iShift+1, jShift+1) * m(iShift+2, jShift+2) * m(iShift+3, jShift+3)
150 - m(iShift+1, jShift+1) * m(iShift+3, jShift+2) * m(iShift+3, jShift+2)
151 - m(iShift+2, jShift+1) * m(iShift+2, jShift+1) * m(iShift+3, jShift+3)
152 + m(iShift+2, jShift+1) * m(iShift+3, jShift+1) * m(iShift+3, jShift+2)
153 + m(iShift+3, jShift+1) * m(iShift+2, jShift+1) * m(iShift+3, jShift+2)
154 - m(iShift+3, jShift+1) * m(iShift+3, jShift+1) * m(iShift+2, jShift+2);
155 // cofactor (1,0) [1]
156 inv(1,0) = -m(iShift+1, jShift+0) * m(iShift+2, jShift+2) * m(iShift+3, jShift+3)
157 + m(iShift+1, jShift+0) * m(iShift+3, jShift+2) * m(iShift+3, jShift+2)
158 + m(iShift+2, jShift+1) * m(iShift+2, jShift+0) * m(iShift+3, jShift+3)
159 - m(iShift+2, jShift+1) * m(iShift+3, jShift+0) * m(iShift+3, jShift+2)
160 - m(iShift+3, jShift+1) * m(iShift+2, jShift+0) * m(iShift+3, jShift+2)
161 + m(iShift+3, jShift+1) * m(iShift+3, jShift+0) * m(iShift+2, jShift+2);
162 // cofactor (2,0) [2]
163 inv(2,0) = m(iShift+1, jShift+0) * m(iShift+2, jShift+1) * m(iShift+3, jShift+3)
164 - m(iShift+1, jShift+0) * m(iShift+3, jShift+1) * m(iShift+3, jShift+2)
165 - m(iShift+1, jShift+1) * m(iShift+2, jShift+0) * m(iShift+3, jShift+3)
166 + m(iShift+1, jShift+1) * m(iShift+3, jShift+0) * m(iShift+3, jShift+2)
167 + m(iShift+3, jShift+1) * m(iShift+2, jShift+0) * m(iShift+3, jShift+1)
168 - m(iShift+3, jShift+1) * m(iShift+3, jShift+0) * m(iShift+2, jShift+1);
169 // cofactor (3,0) [3]
170 inv(3,0) = -m(iShift+1, jShift+0) * m(iShift+2, jShift+1) * m(iShift+3, jShift+2)
171 + m(iShift+1, jShift+0) * m(iShift+3, jShift+1) * m(iShift+2, jShift+2)
172 + m(iShift+1, jShift+1) * m(iShift+2, jShift+0) * m(iShift+3, jShift+2)
173 - m(iShift+1, jShift+1) * m(iShift+3, jShift+0) * m(iShift+2, jShift+2)
174 - m(iShift+2, jShift+1) * m(iShift+2, jShift+0) * m(iShift+3, jShift+1)
175 + m(iShift+2, jShift+1) * m(iShift+3, jShift+0) * m(iShift+2, jShift+1);
176 // cofactor (1,1) [5]
177 inv(1,1) = m(iShift+0, jShift+0) * m(iShift+2, jShift+2) * m(iShift+3, jShift+3)
178 - m(iShift+0, jShift+0) * m(iShift+3, jShift+2) * m(iShift+3, jShift+2)
179 - m(iShift+2, jShift+0) * m(iShift+2, jShift+0) * m(iShift+3, jShift+3)
180 + m(iShift+2, jShift+0) * m(iShift+3, jShift+0) * m(iShift+3, jShift+2)
181 + m(iShift+3, jShift+0) * m(iShift+2, jShift+0) * m(iShift+3, jShift+2)
182 - m(iShift+3, jShift+0) * m(iShift+3, jShift+0) * m(iShift+2, jShift+2);
183 // cofactor (2,1) [6]
184 inv(2,1) = -m(iShift+0, jShift+0) * m(iShift+2, jShift+1) * m(iShift+3, jShift+3)
185 + m(iShift+0, jShift+0) * m(iShift+3, jShift+1) * m(iShift+3, jShift+2)
186 + m(iShift+1, jShift+0) * m(iShift+2, jShift+0) * m(iShift+3, jShift+3)
187 - m(iShift+1, jShift+0) * m(iShift+3, jShift+0) * m(iShift+3, jShift+2)
188 - m(iShift+3, jShift+0) * m(iShift+2, jShift+0) * m(iShift+3, jShift+1)
189 + m(iShift+3, jShift+0) * m(iShift+3, jShift+0) * m(iShift+2, jShift+1);
190 // cofactor (3,1) [7]
191 inv(3,1) = m(iShift+0, jShift+0) * m(iShift+2, jShift+1) * m(iShift+3, jShift+2)
192 - m(iShift+0, jShift+0) * m(iShift+3, jShift+1) * m(iShift+2, jShift+2)
193 - m(iShift+1, jShift+0) * m(iShift+2, jShift+0) * m(iShift+3, jShift+2)
194 + m(iShift+1, jShift+0) * m(iShift+3, jShift+0) * m(iShift+2, jShift+2)
195 + m(iShift+2, jShift+0) * m(iShift+2, jShift+0) * m(iShift+3, jShift+1)
196 - m(iShift+2, jShift+0) * m(iShift+3, jShift+0) * m(iShift+2, jShift+1);
197 // cofactor (2,2) [10]
198 inv(2,2) = m(iShift+0, jShift+0) * m(iShift+1, jShift+1) * m(iShift+3, jShift+3)
199 - m(iShift+0, jShift+0) * m(iShift+3, jShift+1) * m(iShift+3, jShift+1)
200 - m(iShift+1, jShift+0) * m(iShift+1, jShift+0) * m(iShift+3, jShift+3)
201 + m(iShift+1, jShift+0) * m(iShift+3, jShift+0) * m(iShift+3, jShift+1)
202 + m(iShift+3, jShift+0) * m(iShift+1, jShift+0) * m(iShift+3, jShift+1)
203 - m(iShift+3, jShift+0) * m(iShift+3, jShift+0) * m(iShift+1, jShift+1);
204 // cofactor (3,2) [11]
205 inv(3,2) = -m(iShift+0, jShift+0) * m(iShift+1, jShift+1) * m(iShift+3, jShift+2)
206 + m(iShift+0, jShift+0) * m(iShift+3, jShift+1) * m(iShift+2, jShift+1)
207 + m(iShift+1, jShift+0) * m(iShift+1, jShift+0) * m(iShift+3, jShift+2)
208 - m(iShift+1, jShift+0) * m(iShift+3, jShift+0) * m(iShift+2, jShift+1)
209 - m(iShift+2, jShift+0) * m(iShift+1, jShift+0) * m(iShift+3, jShift+1)
210 + m(iShift+2, jShift+0) * m(iShift+3, jShift+0) * m(iShift+1, jShift+1);
211 // cofactor (3,3) [15]
212 inv(3,3) = m(iShift+0, jShift+0) * m(iShift+1, jShift+1) * m(iShift+2, jShift+2)
213 - m(iShift+0, jShift+0) * m(iShift+2, jShift+1) * m(iShift+2, jShift+1)
214 - m(iShift+1, jShift+0) * m(iShift+1, jShift+0) * m(iShift+2, jShift+2)
215 + m(iShift+1, jShift+0) * m(iShift+2, jShift+0) * m(iShift+2, jShift+1)
216 + m(iShift+2, jShift+0) * m(iShift+1, jShift+0) * m(iShift+2, jShift+1)
217 - m(iShift+2, jShift+0) * m(iShift+2, jShift+0) * m(iShift+1, jShift+1);
218 // symmetric
219 inv(0,1) = inv(1,0); inv(0,2) = inv(2,0); inv(0,3) = inv(3,0);
220 inv(1,2) = inv(2,1); inv(1,3) = inv(3,1);
221 inv(2,3) = inv(3,2);
222 // compute determinant
223 Type det = m(iShift+0, jShift+0) * inv(0,0)
224 + m(iShift+1, jShift+0) * inv(1,0)
225 + m(iShift+2, jShift+0) * inv(2,0)
226 + m(iShift+3, jShift+0) * inv(3,0);
227 if (det == Type(0)) return Type(0);
228 // inverse matrix
229 inv /= det;
230 return det;
231 }
239 static Type invertSymMatrixXX( Matrix const& m, CArraySquare<Type, Size_>& inv)
240 {
241#ifdef STKUSELAPACK
242 lapack::SymEigen<Matrix> decomp(m);
243#else
244 SymEigen<Matrix> decomp(m);
245#endif
246 if (!decomp.run()) return Type(0);
247 // compute (generalized) inverse matrix
248 decomp.ginv(inv);
249 return decomp.det();
250 }
251
252};
253
257template<class Matrix, int Size_>
258struct InvertSymMatrixDispatcher
259{
260 typedef typename Matrix::Type Type;
261 inline static Type run( Matrix const& m, CArraySquare<Type, Size_>& inv)
262 {
263 switch (m.sizeRows())
264 {
265 case 1: return InvertSymImpl<Matrix, Size_>::invertSymMatrix11(m, inv);
266 case 2: return InvertSymImpl<Matrix, Size_>::invertSymMatrix22(m, inv);
267 case 3: return InvertSymImpl<Matrix, Size_>::invertSymMatrix33(m, inv);
268 case 4: return InvertSymImpl<Matrix, Size_>::invertSymMatrix44(m, inv);
269 default:return InvertSymImpl<Matrix, Size_>::invertSymMatrixXX(m, inv);
270 }
271 }
272};
273
274template<class Matrix>
275struct InvertSymMatrixDispatcher<Matrix, 1>
276{
277 typedef typename Matrix::Type Type;
278 inline static Type run( Matrix const& m, CArraySquare<Type, 1>& inv)
279 { return InvertSymImpl<Matrix, 1>::invertSymMatrix11(m, inv);}
280};
281
282template<class Matrix>
283struct InvertSymMatrixDispatcher<Matrix, 2>
284{
285 typedef typename Matrix::Type Type;
286 inline static Type run( Matrix const& m, CArraySquare<Type, 2>& inv)
287 { return InvertSymImpl<Matrix, 2>::invertSymMatrix22(m, inv);}
288};
289
290template<class Matrix>
291struct InvertSymMatrixDispatcher<Matrix, 3>
292{
293 typedef typename Matrix::Type Type;
294 inline static Type run( Matrix const& m, CArraySquare<Type, 3>& inv)
295 { return InvertSymImpl<Matrix, 3>::invertSymMatrix33(m, inv);}
296};
297
298template<class Matrix>
299struct InvertSymMatrixDispatcher<Matrix, 4>
300{
301 typedef typename Matrix::Type Type;
302 inline static Type run( Matrix const& m, CArraySquare<Type, 4>& inv)
303 { return InvertSymImpl<Matrix, 4>::invertSymMatrix44(m, inv);}
304};
305
306} // namespace hidden
307
308} // namespace STK
309
310#endif /* STK_SYMINVERT_H */
STK_SymEigen.h
In this file we define the SymEigen class (for a symmetric matrix).
STK_lapack_SymEigen.h
In this file we define the enclosing class of the syevr lapck routine.
STK::ISymEigen::ginv
ArraySquare & ginv(ArraySquare &res) const
Compute the generalized inverse of the symmetric matrix and put the result in res.
Definition STK_ISymEigen.h:156
STK::ISymEigen::run
virtual bool run()
Find the eigenvalues and eigenvectors of the matrix.
Definition STK_ISymEigen.h:189
STK::ISymEigen::det
Type const & det() const
Definition STK_ISymEigen.h:127
STK::MultidimRegression
The MultidimRegression class allows to regress a multidimensional output variable among a multivariat...
Definition STK_MultidimRegression.h:52
STK::hidden::InvertSymImpl::invertSymMatrix33
static Type invertSymMatrix33(Matrix const &m, CArraySquare< Type, Size_ > &inv)
compute the inverse of the symmetric 3x3 matrix m using only its lower part and store the result in i...
Definition STK_InvertSym.h:108
STK::hidden::InvertSymImpl::invertSymMatrix44
static Type invertSymMatrix44(Matrix const &m, CArraySquare< Type, Size_ > &inv)
compute the inverse of the symmetric 4x4 matrix m using only its lower part and store the result in i...
Definition STK_InvertSym.h:144
STK::hidden::InvertSymImpl::invertSymMatrixXX
static Type invertSymMatrixXX(Matrix const &m, CArraySquare< Type, Size_ > &inv)
compute the inverse of the symmetric 4x4 matrix m using only its lower part and store the result in i...
Definition STK_InvertSym.h:239
STK::hidden::InvertSymImpl::invertSymMatrix11
static Type invertSymMatrix11(Matrix const &m, CArraySquare< Type, Size_ > &inv)
compute the inverse of the symmetric matrix m of size 1x1 and store the result in inv.
Definition STK_InvertSym.h:62
STK::hidden::InvertSymImpl::invertSymMatrix22
static Type invertSymMatrix22(Matrix const &m, CArraySquare< Type, Size_ > &inv)
compute the inverse of the symmetric 2x2 matrix m using only its lower part and store the result in i...
Definition STK_InvertSym.h:82
STK
The namespace STK is the main domain space of the Statistical ToolKit project.
STK::hidden::InvertSymImpl
Implementation of the computation of the inverse of a lower triangular matrix.
Definition STK_InvertSym.h:54
STK::hidden::InvertSymMatrixDispatcher< Matrix, 1 >::run
static Type run(Matrix const &m, CArraySquare< Type, 1 > &inv)
Definition STK_InvertSym.h:278
STK::hidden::InvertSymMatrixDispatcher< Matrix, 2 >::run
static Type run(Matrix const &m, CArraySquare< Type, 2 > &inv)
Definition STK_InvertSym.h:286
STK::hidden::InvertSymMatrixDispatcher< Matrix, 3 >::run
static Type run(Matrix const &m, CArraySquare< Type, 3 > &inv)
Definition STK_InvertSym.h:294
STK::hidden::InvertSymMatrixDispatcher< Matrix, 4 >::run
static Type run(Matrix const &m, CArraySquare< Type, 4 > &inv)
Definition STK_InvertSym.h:302
STK::hidden::InvertSymMatrixDispatcher
computing the inverse of a symmetric matrix.
Definition STK_InvertSym.h:259
STK::hidden::InvertSymMatrixDispatcher::run
static Type run(Matrix const &m, CArraySquare< Type, Size_ > &inv)
Definition STK_InvertSym.h:261