solveSystem.h
Go to the documentation of this file.
1 //===========================================================================
2 /*!
3  *
4  *
5  * \brief Some operations for matrices.
6  *
7  *
8  *
9  *
10  * \author O. Krause
11  * \date 2011
12  *
13  *
14  * \par Copyright 1995-2015 Shark Development Team
15  *
16  * <BR><HR>
17  * This file is part of Shark.
18  * <http://image.diku.dk/shark/>
19  *
20  * Shark is free software: you can redistribute it and/or modify
21  * it under the terms of the GNU Lesser General Public License as published
22  * by the Free Software Foundation, either version 3 of the License, or
23  * (at your option) any later version.
24  *
25  * Shark is distributed in the hope that it will be useful,
26  * but WITHOUT ANY WARRANTY; without even the implied warranty of
27  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
28  * GNU Lesser General Public License for more details.
29  *
30  * You should have received a copy of the GNU Lesser General Public License
31  * along with Shark. If not, see <http://www.gnu.org/licenses/>.
32  *
33  */
34 //===========================================================================
35 
36 #ifndef SHARK_LINALG_SOLVE_SYSTEM_H
37 #define SHARK_LINALG_SOLVE_SYSTEM_H
38 
40 
41 namespace shark{ namespace blas{
42 
43 /**
44  * \ingroup shark_globals
45  *
46  * @{
47  */
48 
49 /// \brief In-Place System of linear equations solver.
50 ///
51 ///Solves a system of linear equations
52 ///Ax=b
53 ///for x, using LU decomposition and
54 ///backward substitution sotring the results in b.
55 ///This Method is in
56 ///no way optimized for sparse matrices.
57 ///Be aware, that the matrix must have full rank!
58 template<class MatT,class VecT>
60  matrix_expression<MatT> const& A,
61  vector_expression<VecT>& b
62 );
63 /// \brief System of linear equations solver.
64 ///
65 /// Solves asystem of linear equations
66 /// Ax=b
67 /// for x, using LU decomposition and
68 /// backward substitution. This Method is in
69 /// no way optimized for sparse matrices.
70 /// Be aware, that the matrix must have full rank!
71 template<class MatT,class Vec1T,class Vec2T>
72 void solveSystem(
73  const matrix_expression<MatT> & A,
74  vector_expression<Vec1T>& x,
75  const vector_expression<Vec2T> & b
76 );
77 
78 /// \brief In-Place system of linear equations solver.
79 ///
80 ///Solves multiple systems of linear equations
81 ///Ax_1=b_1
82 ///Ax_2=b_2
83 ///...
84 ///=>AX=B
85 ///for X, using LU decomposition and
86 ///backward substitution and stores the result in b
87 ///Note, that B=(b_1,...,b_n), so the right hand sides are stored as columns
88 ///This Method is in no way optimized for sparse matrices.
89 ///Be aware, that the matrix must have full rank!
90 template<class MatT,class Mat2T>
92  matrix_expression<MatT> const& A,
93  matrix_expression<Mat2T>& B
94 );
95 /// \brief System of linear equations solver.
96 ///
97 /// Solves multiple systems of linear equations
98 /// Ax_1=b_1
99 /// Ax_2=b_2
100 /// ...
101 /// =>AX=B
102 /// for X, using LU decomposition and
103 /// backward substitution.
104 /// Note, that B=(b_1,...,b_n), so the right hand sides are stored as columns
105 /// This Method is in no way optimized for sparse matrices.
106 /// Be aware that the matrix must have full rank!
107 template<class MatT,class Mat1T,class Mat2T>
108 void solveSystem(
109  const matrix_expression<MatT> & A,
110  matrix_expression<Mat1T>& X,
111  const matrix_expression<Mat2T> & B
112 );
113 
114 /// \brief System of symmetric linear equations solver. The result is stored in b
115 ///
116 /// Solves a system of linear equations
117 /// Ax=b
118 /// for x, using Cholesky decomposition and
119 /// backward substitution. and stores the result in b.
120 /// A must be symmetric.
121 /// This method is in no way optimized for sparse matrices.
122 /// Be aware, that the matrix must have full rank!
123 template<class System, class MatT,class VecT>
125  matrix_expression<MatT> const& A,
126  vector_expression<VecT>& b
127 );
128 
129 /// \brief System of symmetric linear equations solver.
130 ///
131 /// Solves multiple systems of linear equations
132 /// Ax_1=b_1
133 /// Ax_1=b_2
134 /// ...
135 /// =>AX=B
136 /// or XA=B
137 /// for X, using cholesky decomposition and
138 /// backward substitution. The first template parameter is used
139 /// to decide which type of system is solved
140 /// Note, that B=(b_1,...,b_n), so the right hand sides are stored as columns.
141 /// A must be symmetric.
142 /// This Method is in no way optimized for sparse matrices.
143 /// Be aware, that the matrix must have full rank!
144 /// Also the result is stored in B directly so it"s contents are destroyed.
145 /// @param A the system matrix A
146 /// @param B the right hand side of the LGS, also stores the result
147 template<class System, class MatT,class Mat1T>
149  matrix_expression<MatT> const& A,
150  matrix_expression<Mat1T>& B
151 );
152 
153 /// \brief System of symmetric linear equations solver.
154 ///
155 /// Solves a system of linear equations
156 /// Ax=b
157 /// for x, using Cholesky decomposition and
158 /// backward substitution. A must be symmetric.
159 /// This Method is in no way optimized for sparse matrices.
160 /// Be aware, that the matrix must have full rank!
161 template<class System,class MatT,class Vec1T,class Vec2T>
163  matrix_expression<MatT> const& A,
164  vector_expression<Vec1T>& x,
165  vector_expression<Vec2T> const& b
166 );
167 /// \brief System of symmetric linear equations solver.
168 ///
169 /// Solves multiple systems of linear equations
170 /// Ax_1=b_1
171 /// Ax_1=b_2
172 /// ...
173 /// =>AX=B
174 /// or XA = B
175 /// for X, using cholesky decomposition and
176 /// backward substitution. The first template parameter is used
177 /// to decide which type of system is solved
178 /// Note, that B=(b_1,...,b_n), so the right hand sides are stored as columns.
179 /// A must be symmetric.
180 /// This Method is in no way optimized for sparse matrices.
181 /// Be aware, that the matrix must have full rank!
182 /// @param A the system matrix A
183 /// @param X the stored result of the solution of LGS
184 /// @param B the right hand side of the LGS
185 template<class System,class MatT,class Mat1T,class Mat2T>
187  matrix_expression<MatT> const& A,
188  matrix_expression<Mat1T>& X,
189  matrix_expression<Mat2T> const& B
190 );
191 
192 
193 /// \brief Solves a square system of linear equations without full rank.
194 ///
195 /// Solves the system Ax= b or x^TA=b^T when A is
196 /// symmetric positive semi-definite.
197 /// If b is not in the span of Ax or xA, the least squares solution is used,
198 /// that is we minimize ||Ax-b||^2
199 ///
200 /// The computation is carried out in-place.
201 /// The algorithm can be looked up in
202 /// "Fast Computation of Moore-Penrose Inverse Matrices"
203 /// Pierre Courrieu, 2005
204 ///
205 /// \param A \f$ n \times n \f$ input matrix.
206 /// \param b right hand side vector.
207 template<class System,class MatT,class VecT>
209  matrix_expression<MatT> const& A,
210  vector_expression<VecT>& b
211 );
212 
213 /// \brief Solves multiple square system of linear equations without full rank.
214 ///
215 /// Solves multiple systems of linear equations
216 /// Ax_1=b_1
217 /// Ax_1=b_2
218 /// ...
219 /// =>AX=B
220 /// or XA = B
221 /// A must be symmetric positive semi-definite - thus is not required to have full rank.
222 /// Note, that B=(b_1,...,b_n), so the right hand sides are stored as columns.
223 /// If the b_i are not in the span of Ax_i or x_i^TA, the least squares solution is used,
224 /// that is we minimize ||Ax_i-b_i||^2
225 ///
226 /// The computation is carried out in-place.
227 /// The algorithm can be looked up in
228 /// "Fast Computation of Moore-Penrose Inverse Matrices"
229 /// Pierre Courrieu, 2005
230 ///
231 /// \param A \f$ n \times n \f$ input matrix.
232 /// \param B \f$ n \times k \f$ right hand side matrix.
233 template<class System,class Mat1T,class Mat2T>
235  matrix_expression<Mat1T> const& A,
236  matrix_expression<Mat2T>& B
237 );
238 
239 /// \brief Solves a non-square system of linear equations.
240 ///
241 /// Given a \f$ m \times n \f$ input matrix A this function uses
242 /// the generalized inverse of A to solve the system of linear equations.
243 /// If b is not in the span of Ax or xA, the least squares solution is used,
244 /// that is we minimize ||Ax-b||^2
245 ///
246 /// The computation is carried out in-place.
247 ///
248 /// \param A \f$ n \times m \f$ input matrix.
249 /// \param b right hand side of the problem.
250 template<class System,class MatT,class VecT>
252  matrix_expression<MatT> const& A,
253  vector_expression<VecT>& b
254 );
255 
256 
257 /// \brief Solves multiple non-square systems of linear equations.
258 ///
259 /// Given a \f$ m \times n \f$ input matrix A this function uses
260 /// the generalized inverse of A to solve the system of linear equations AX=B.
261 /// If b_i is not in the span of Ax_i or x_iA, the least squares solution is used,
262 /// that is we minimize ||Ax_i-b_i||^2 for all columns b_i of B.
263 ///
264 /// The computation is carried out in-place.
265 ///
266 /// \param A \f$ n \times m \f$ input matrix.
267 /// \param B \f$ n \times k \f$ right hand sied matrix.
268 template<class System,class MatA,class MatB>
270  matrix_expression<MatA> const& A,
271  matrix_expression<MatB>& B
272 );
273 
274 /// \brief Approximates the solution of a linear system of equation Ax=b.
275 ///
276 /// Most often there is no need for the exact solution of a system of linear
277 /// equations. Instead only a good approximation needs to be found.
278 /// In this case an iterative method can be used which stops when
279 /// a suitable exact solution is found. For a lot of systems this already happens
280 /// after a very low number of iterations.
281 /// Every iteration has complexity O(n^2) and after n iterations the
282 /// exact solution is found. However if this solution is needed, the other
283 /// methods, as for example solveSymmPosDefSystem are more suitable.
284 ///
285 /// This algorithm does not require A to have full rank, however it must be
286 /// positive semi-definite.
287 ///
288 /// This algorithm stops after the maximum number of iterations is
289 /// exceeded or after the max-norm of the residual \f$ r_k= -Ax_k+b\f$ is
290 /// smaller than epsilon.
291 ///
292 /// The initial solution argument governs whether x already stores a possible starting point.
293 /// If this is true, it is checked whether it is better than
294 /// starting from 0 (i.e. the max-norm of the initial residual is smaller than -b).
295 ///
296 /// \param A the positive semi-definite n x n-Matrix
297 /// \param x the solution vector
298 /// \param b the right hand side
299 /// \param epsilon stopping criterium for the residual
300 /// \param maxIterations the maximum number of iterations
301 /// \param initialSolution if this is true, x stores an initial guess of the solution
302 template<class MatT, class VecT, class VecT2>
304  matrix_expression<MatT> const& A,
306  vector_expression<VecT2> const& b,
307  double epsilon = 1.e-10,
308  bool initialSolution = false,
309  unsigned int maxIterations = 0
310 ){
311  SIZE_CHECK(A().size1()==A().size2());
312  SIZE_CHECK(A().size1()==b().size());
313 
314  std::size_t dim = b().size();
315  std::size_t maxIt = (maxIterations == 0)? dim: maxIterations;
316 
317  typedef typename VecT::value_type value_type;
318  vector<value_type> r = b;//current residual
319  if(initialSolution){
320  SIZE_CHECK(x().size() == dim);
321  noalias(r) -= prod(A,x);
322  if(norm_inf(r) > norm_inf(b)){
323  x().clear();
324  r = b;
325  }
326  }
327  else{
328  ensure_size(x,dim);
329  x().clear();
330  }
331 
332  vector<value_type> rnext(dim); //the next residual
333  vector<value_type> p = r; //the search direction- initially it is the gradient direction
334  vector<value_type> Ap(dim); //stores prod(A,p)
335 
336  for(std::size_t i = 0; i != maxIt; ++i){
337  noalias(Ap) = prod(A,p);
338  double rsqr=inner_prod(r,r);
339  double alpha = rsqr/inner_prod(p,Ap);
340  noalias(x())+=alpha*p;
341  noalias(rnext) = r - alpha * Ap;
342  if(norm_inf(rnext)<epsilon)
343  break;
344 
345  double beta = inner_prod(rnext,rnext)/rsqr;
346  p*=beta;
347  noalias(p) +=rnext;
348  swap(r,rnext);
349  }
350 }
351 
352 /// \brief Approximates the solution of a linear system of equation Ax=b, storing the solution in b.
353 ///
354 /// Most often there is no need for the exact solution of a system of linear
355 /// equations. Instead only a good approximation needs to be found.
356 /// In this case an iterative method can be used which stops when
357 /// a suitable exact solution is found. For a lot of systems this already happens
358 /// after a very low number of iterations.
359 /// Every iteration has complexity O(n^2) and after n iterations the
360 /// exact solution is found. However if this solution is needed, the other
361 /// methods, as for xample solveSymmPosDefSystem are more suitable.
362 ///
363 /// This algorithm stops after the maximum number of iterations is
364 /// exceeded or after the max-norm of the residual \f$ r_k= Ax_k-b\f$ is
365 /// smaller than epsilon. The reuslt is stored in b afterwars
366 ///
367 /// \param A the positive semi-definite n x n-Matrix
368 /// \param b the right hand side which also stores the final solution
369 /// \param epsilon stopping criterium for the residual
370 /// \param maxIterations the maximum number of iterations
371 template<class MatT, class VecT>
373  matrix_expression<MatT> const& A,
375  double epsilon = 1.e-10,
376  unsigned int maxIterations = 0
377 ){
378  SIZE_CHECK(A().size1()==A().size2());
379  SIZE_CHECK(A().size1()==b().size());
380  vector<typename VecT::value_type> x(b.size(),0.0);
381  approxsolveSymmPosDefSystem(A,x,b,epsilon,false,maxIterations);
382  swap(x,b);
383 }
384 
385 /** @}*/
386 }}
387 
388 
389 #include "Impl/solveSystem.inl"
390 #endif