[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]

linear_solve.hxx
1 /************************************************************************/
2 /* */
3 /* Copyright 2003-2008 by Gunnar Kedenburg and Ullrich Koethe */
4 /* */
5 /* This file is part of the VIGRA computer vision library. */
6 /* The VIGRA Website is */
7 /* http://hci.iwr.uni-heidelberg.de/vigra/ */
8 /* Please direct questions, bug reports, and contributions to */
9 /* ullrich.koethe@iwr.uni-heidelberg.de or */
10 /* vigra@informatik.uni-hamburg.de */
11 /* */
12 /* Permission is hereby granted, free of charge, to any person */
13 /* obtaining a copy of this software and associated documentation */
14 /* files (the "Software"), to deal in the Software without */
15 /* restriction, including without limitation the rights to use, */
16 /* copy, modify, merge, publish, distribute, sublicense, and/or */
17 /* sell copies of the Software, and to permit persons to whom the */
18 /* Software is furnished to do so, subject to the following */
19 /* conditions: */
20 /* */
21 /* The above copyright notice and this permission notice shall be */
22 /* included in all copies or substantial portions of the */
23 /* Software. */
24 /* */
25 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */
26 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */
27 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */
28 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */
29 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */
30 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */
31 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */
32 /* OTHER DEALINGS IN THE SOFTWARE. */
33 /* */
34 /************************************************************************/
35 
36 
37 #ifndef VIGRA_LINEAR_SOLVE_HXX
38 #define VIGRA_LINEAR_SOLVE_HXX
39 
40 #include <ctype.h>
41 #include <string>
42 #include "mathutil.hxx"
43 #include "matrix.hxx"
44 #include "singular_value_decomposition.hxx"
45 
46 
47 namespace vigra
48 {
49 
50 namespace linalg
51 {
52 
53 namespace detail {
54 
55 template <class T, class C1>
56 T determinantByLUDecomposition(MultiArrayView<2, T, C1> const & a)
57 {
58  typedef MultiArrayShape<2>::type Shape;
59 
60  MultiArrayIndex m = rowCount(a), n = columnCount(a);
61  vigra_precondition(n == m,
62  "determinant(): square matrix required.");
63 
64  Matrix<T> LU(a);
65  T det = 1.0;
66 
67  for (MultiArrayIndex j = 0; j < n; ++j)
68  {
69  // Apply previous transformations.
70  for (MultiArrayIndex i = 0; i < m; ++i)
71  {
72  MultiArrayIndex end = std::min(i, j);
73  T s = dot(rowVector(LU, Shape(i,0), end), columnVector(LU, Shape(0,j), end));
74  LU(i,j) = LU(i,j) -= s;
75  }
76 
77  // Find pivot and exchange if necessary.
78  MultiArrayIndex p = j + argMax(abs(columnVector(LU, Shape(j,j), m)));
79  if (p != j)
80  {
81  rowVector(LU, p).swapData(rowVector(LU, j));
82  det = -det;
83  }
84 
85  det *= LU(j,j);
86 
87  // Compute multipliers.
88  if (LU(j,j) != 0.0)
89  columnVector(LU, Shape(j+1,j), m) /= LU(j,j);
90  else
91  break; // det is zero
92  }
93  return det;
94 }
95 
96 // returns the new value of 'a' (when this Givens rotation is applied to 'a' and 'b')
97 // the new value of 'b' is zero, of course
98 template <class T>
99 T givensCoefficients(T a, T b, T & c, T & s)
100 {
101  if(abs(a) < abs(b))
102  {
103  T t = a/b,
104  r = std::sqrt(1.0 + t*t);
105  s = 1.0 / r;
106  c = t*s;
107  return r*b;
108  }
109  else if(a != 0.0)
110  {
111  T t = b/a,
112  r = std::sqrt(1.0 + t*t);
113  c = 1.0 / r;
114  s = t*c;
115  return r*a;
116  }
117  else // a == b == 0.0
118  {
119  c = 1.0;
120  s = 0.0;
121  return 0.0;
122  }
123 }
124 
125 // see Golub, van Loan: Algorithm 5.1.3 (p. 216)
126 template <class T>
127 bool givensRotationMatrix(T a, T b, Matrix<T> & gTranspose)
128 {
129  if(b == 0.0)
130  return false; // no rotation needed
131  givensCoefficients(a, b, gTranspose(0,0), gTranspose(0,1));
132  gTranspose(1,1) = gTranspose(0,0);
133  gTranspose(1,0) = -gTranspose(0,1);
134  return true;
135 }
136 
137 // reflections are symmetric matrices and can thus be applied to rows
138 // and columns in the same way => code simplification relative to rotations
139 template <class T>
140 inline bool
141 givensReflectionMatrix(T a, T b, Matrix<T> & g)
142 {
143  if(b == 0.0)
144  return false; // no reflection needed
145  givensCoefficients(a, b, g(0,0), g(0,1));
146  g(1,1) = -g(0,0);
147  g(1,0) = g(0,1);
148  return true;
149 }
150 
151 // see Golub, van Loan: Algorithm 5.2.2 (p. 227) and Section 12.5.2 (p. 608)
152 template <class T, class C1, class C2>
153 bool
154 qrGivensStepImpl(MultiArrayIndex i, MultiArrayView<2, T, C1> r, MultiArrayView<2, T, C2> rhs)
155 {
156  typedef typename Matrix<T>::difference_type Shape;
157 
158  const MultiArrayIndex m = rowCount(r);
159  const MultiArrayIndex n = columnCount(r);
160  const MultiArrayIndex rhsCount = columnCount(rhs);
161  vigra_precondition(m == rowCount(rhs),
162  "qrGivensStepImpl(): Matrix shape mismatch.");
163 
164  Matrix<T> givens(2,2);
165  for(int k=m-1; k>(int)i; --k)
166  {
167  if(!givensReflectionMatrix(r(k-1,i), r(k,i), givens))
168  continue; // r(k,i) was already zero
169 
170  r(k-1,i) = givens(0,0)*r(k-1,i) + givens(0,1)*r(k,i);
171  r(k,i) = 0.0;
172 
173  r.subarray(Shape(k-1,i+1), Shape(k+1,n)) = givens*r.subarray(Shape(k-1,i+1), Shape(k+1,n));
174  rhs.subarray(Shape(k-1,0), Shape(k+1,rhsCount)) = givens*rhs.subarray(Shape(k-1,0), Shape(k+1,rhsCount));
175  }
176  return r(i,i) != 0.0;
177 }
178 
179 // see Golub, van Loan: Section 12.5.2 (p. 608)
180 template <class T, class C1, class C2, class Permutation>
181 void
182 upperTriangularCyclicShiftColumns(MultiArrayIndex i, MultiArrayIndex j,
183  MultiArrayView<2, T, C1> &r, MultiArrayView<2, T, C2> &rhs, Permutation & permutation)
184 {
185  typedef typename Matrix<T>::difference_type Shape;
186 
187  const MultiArrayIndex m = rowCount(r);
188  const MultiArrayIndex n = columnCount(r);
189  const MultiArrayIndex rhsCount = columnCount(rhs);
190  vigra_precondition(i < n && j < n,
191  "upperTriangularCyclicShiftColumns(): Shift indices out of range.");
192  vigra_precondition(m == rowCount(rhs),
193  "upperTriangularCyclicShiftColumns(): Matrix shape mismatch.");
194 
195  if(j == i)
196  return;
197  if(j < i)
198  std::swap(j,i);
199 
200  Matrix<T> t = columnVector(r, i);
201  MultiArrayIndex ti = permutation[i];
202  for(MultiArrayIndex k=i; k<j;++k)
203  {
204  columnVector(r, k) = columnVector(r, k+1);
205  permutation[k] = permutation[k+1];
206  }
207  columnVector(r, j) = t;
208  permutation[j] = ti;
209 
210  Matrix<T> givens(2,2);
211  for(MultiArrayIndex k=i; k<j; ++k)
212  {
213  if(!givensReflectionMatrix(r(k,k), r(k+1,k), givens))
214  continue; // r(k+1,k) was already zero
215 
216  r(k,k) = givens(0,0)*r(k,k) + givens(0,1)*r(k+1,k);
217  r(k+1,k) = 0.0;
218 
219  r.subarray(Shape(k,k+1), Shape(k+2,n)) = givens*r.subarray(Shape(k,k+1), Shape(k+2,n));
220  rhs.subarray(Shape(k,0), Shape(k+2,rhsCount)) = givens*rhs.subarray(Shape(k,0), Shape(k+2,rhsCount));
221  }
222 }
223 
224 // see Golub, van Loan: Section 12.5.2 (p. 608)
225 template <class T, class C1, class C2, class Permutation>
226 void
227 upperTriangularSwapColumns(MultiArrayIndex i, MultiArrayIndex j,
228  MultiArrayView<2, T, C1> &r, MultiArrayView<2, T, C2> &rhs, Permutation & permutation)
229 {
230  typedef typename Matrix<T>::difference_type Shape;
231 
232  const MultiArrayIndex m = rowCount(r);
233  const MultiArrayIndex n = columnCount(r);
234  const MultiArrayIndex rhsCount = columnCount(rhs);
235  vigra_precondition(i < n && j < n,
236  "upperTriangularSwapColumns(): Swap indices out of range.");
237  vigra_precondition(m == rowCount(rhs),
238  "upperTriangularSwapColumns(): Matrix shape mismatch.");
239 
240  if(j == i)
241  return;
242  if(j < i)
243  std::swap(j,i);
244 
245  columnVector(r, i).swapData(columnVector(r, j));
246  std::swap(permutation[i], permutation[j]);
247 
248  Matrix<T> givens(2,2);
249  for(int k=m-1; k>(int)i; --k)
250  {
251  if(!givensReflectionMatrix(r(k-1,i), r(k,i), givens))
252  continue; // r(k,i) was already zero
253 
254  r(k-1,i) = givens(0,0)*r(k-1,i) + givens(0,1)*r(k,i);
255  r(k,i) = 0.0;
256 
257  r.subarray(Shape(k-1,i+1), Shape(k+1,n)) = givens*r.subarray(Shape(k-1,i+1), Shape(k+1,n));
258  rhs.subarray(Shape(k-1,0), Shape(k+1,rhsCount)) = givens*rhs.subarray(Shape(k-1,0), Shape(k+1,rhsCount));
259  }
260  MultiArrayIndex end = std::min(j, m-1);
261  for(MultiArrayIndex k=i+1; k<end; ++k)
262  {
263  if(!givensReflectionMatrix(r(k,k), r(k+1,k), givens))
264  continue; // r(k+1,k) was already zero
265 
266  r(k,k) = givens(0,0)*r(k,k) + givens(0,1)*r(k+1,k);
267  r(k+1,k) = 0.0;
268 
269  r.subarray(Shape(k,k+1), Shape(k+2,n)) = givens*r.subarray(Shape(k,k+1), Shape(k+2,n));
270  rhs.subarray(Shape(k,0), Shape(k+2,rhsCount)) = givens*rhs.subarray(Shape(k,0), Shape(k+2,rhsCount));
271  }
272 }
273 
274 // see Lawson & Hanson: Algorithm H1 (p. 57)
275 template <class T, class C1, class C2, class U>
276 bool householderVector(MultiArrayView<2, T, C1> const & v, MultiArrayView<2, T, C2> & u, U & vnorm)
277 {
278  vnorm = (v(0,0) > 0.0)
279  ? -norm(v)
280  : norm(v);
281  U f = std::sqrt(vnorm*(vnorm - v(0,0)));
282 
283  if(f == NumericTraits<U>::zero())
284  {
285  u.init(NumericTraits<T>::zero());
286  return false;
287  }
288  else
289  {
290  u(0,0) = (v(0,0) - vnorm) / f;
291  for(MultiArrayIndex k=1; k<rowCount(u); ++k)
292  u(k,0) = v(k,0) / f;
293  return true;
294  }
295 }
296 
297 // see Lawson & Hanson: Algorithm H1 (p. 57)
298 template <class T, class C1, class C2, class C3>
299 bool
300 qrHouseholderStepImpl(MultiArrayIndex i, MultiArrayView<2, T, C1> & r,
301  MultiArrayView<2, T, C2> & rhs, MultiArrayView<2, T, C3> & householderMatrix)
302 {
303  typedef typename Matrix<T>::difference_type Shape;
304 
305  const MultiArrayIndex m = rowCount(r);
306  const MultiArrayIndex n = columnCount(r);
307  const MultiArrayIndex rhsCount = columnCount(rhs);
308 
309  vigra_precondition(i < n && i < m,
310  "qrHouseholderStepImpl(): Index i out of range.");
311 
312  Matrix<T> u(m-i,1);
313  T vnorm;
314  bool nontrivial = householderVector(columnVector(r, Shape(i,i), m), u, vnorm);
315 
316  r(i,i) = vnorm;
317  columnVector(r, Shape(i+1,i), m).init(NumericTraits<T>::zero());
318 
319  if(columnCount(householderMatrix) == n)
320  columnVector(householderMatrix, Shape(i,i), m) = u;
321 
322  if(nontrivial)
323  {
324  for(MultiArrayIndex k=i+1; k<n; ++k)
325  columnVector(r, Shape(i,k), m) -= dot(columnVector(r, Shape(i,k), m), u) * u;
326  for(MultiArrayIndex k=0; k<rhsCount; ++k)
327  columnVector(rhs, Shape(i,k), m) -= dot(columnVector(rhs, Shape(i,k), m), u) * u;
328  }
329  return r(i,i) != 0.0;
330 }
331 
332 template <class T, class C1, class C2>
333 bool
334 qrColumnHouseholderStep(MultiArrayIndex i, MultiArrayView<2, T, C1> &r, MultiArrayView<2, T, C2> &rhs)
335 {
336  Matrix<T> dontStoreHouseholderVectors; // intentionally empty
337  return qrHouseholderStepImpl(i, r, rhs, dontStoreHouseholderVectors);
338 }
339 
340 template <class T, class C1, class C2>
341 bool
342 qrRowHouseholderStep(MultiArrayIndex i, MultiArrayView<2, T, C1> &r, MultiArrayView<2, T, C2> & householderMatrix)
343 {
344  Matrix<T> dontTransformRHS; // intentionally empty
345  MultiArrayView<2, T, StridedArrayTag> rt = transpose(r),
346  ht = transpose(householderMatrix);
347  return qrHouseholderStepImpl(i, rt, dontTransformRHS, ht);
348 }
349 
350 // O(n) algorithm due to Bischof: Incremental Condition Estimation, 1990
351 template <class T, class C1, class C2, class SNType>
352 void
353 incrementalMaxSingularValueApproximation(MultiArrayView<2, T, C1> const & newColumn,
354  MultiArrayView<2, T, C2> & z, SNType & v)
355 {
356  typedef typename Matrix<T>::difference_type Shape;
357  MultiArrayIndex n = rowCount(newColumn) - 1;
358 
359  SNType vneu = squaredNorm(newColumn);
360  T yv = dot(columnVector(newColumn, Shape(0,0),n), columnVector(z, Shape(0,0),n));
361  // use atan2 as it is robust against overflow/underflow
362  T t = 0.5*std::atan2(T(2.0*yv), T(sq(v)-vneu)),
363  s = std::sin(t),
364  c = std::cos(t);
365  v = std::sqrt(sq(c*v) + sq(s)*vneu + 2.0*s*c*yv);
366  columnVector(z, Shape(0,0),n) = c*columnVector(z, Shape(0,0),n) + s*columnVector(newColumn, Shape(0,0),n);
367  z(n,0) = s*newColumn(n,0);
368 }
369 
370 // O(n) algorithm due to Bischof: Incremental Condition Estimation, 1990
371 template <class T, class C1, class C2, class SNType>
372 void
373 incrementalMinSingularValueApproximation(MultiArrayView<2, T, C1> const & newColumn,
374  MultiArrayView<2, T, C2> & z, SNType & v, double tolerance)
375 {
376  typedef typename Matrix<T>::difference_type Shape;
377 
378  if(v <= tolerance)
379  {
380  v = 0.0;
381  return;
382  }
383 
384  MultiArrayIndex n = rowCount(newColumn) - 1;
385 
386  T gamma = newColumn(n,0);
387  if(gamma == 0.0)
388  {
389  v = 0.0;
390  return;
391  }
392 
393  T yv = dot(columnVector(newColumn, Shape(0,0),n), columnVector(z, Shape(0,0),n));
394  // use atan2 as it is robust against overflow/underflow
395  T t = 0.5*std::atan2(T(-2.0*yv), T(squaredNorm(gamma / v) + squaredNorm(yv) - 1.0)),
396  s = std::sin(t),
397  c = std::cos(t);
398  columnVector(z, Shape(0,0),n) *= c;
399  z(n,0) = (s - c*yv) / gamma;
400  v *= norm(gamma) / hypot(c*gamma, v*(s - c*yv));
401 }
402 
403 // QR algorithm with optional column pivoting
404 template <class T, class C1, class C2, class C3>
405 unsigned int
406 qrTransformToTriangularImpl(MultiArrayView<2, T, C1> & r, MultiArrayView<2, T, C2> & rhs, MultiArrayView<2, T, C3> & householder,
407  ArrayVector<MultiArrayIndex> & permutation, double epsilon)
408 {
409  typedef typename Matrix<T>::difference_type Shape;
410  typedef typename NormTraits<MultiArrayView<2, T, C1> >::NormType NormType;
411  typedef typename NormTraits<MultiArrayView<2, T, C1> >::SquaredNormType SNType;
412 
413  const MultiArrayIndex m = rowCount(r);
414  const MultiArrayIndex n = columnCount(r);
415  const MultiArrayIndex maxRank = std::min(m, n);
416 
417  vigra_precondition(m >= n,
418  "qrTransformToTriangularImpl(): Coefficient matrix with at least as many rows as columns required.");
419 
420  const MultiArrayIndex rhsCount = columnCount(rhs);
421  bool transformRHS = rhsCount > 0;
422  vigra_precondition(!transformRHS || m == rowCount(rhs),
423  "qrTransformToTriangularImpl(): RHS matrix shape mismatch.");
424 
425  bool storeHouseholderSteps = columnCount(householder) > 0;
426  vigra_precondition(!storeHouseholderSteps || r.shape() == householder.shape(),
427  "qrTransformToTriangularImpl(): Householder matrix shape mismatch.");
428 
429  bool pivoting = permutation.size() > 0;
430  vigra_precondition(!pivoting || n == (MultiArrayIndex)permutation.size(),
431  "qrTransformToTriangularImpl(): Permutation array size mismatch.");
432 
433  if(n == 0)
434  return 0; // trivial solution
435 
436  Matrix<SNType> columnSquaredNorms;
437  if(pivoting)
438  {
439  columnSquaredNorms.reshape(Shape(1,n));
440  for(MultiArrayIndex k=0; k<n; ++k)
441  columnSquaredNorms[k] = squaredNorm(columnVector(r, k));
442 
443  int pivot = argMax(columnSquaredNorms);
444  if(pivot != 0)
445  {
446  columnVector(r, 0).swapData(columnVector(r, pivot));
447  std::swap(columnSquaredNorms[0], columnSquaredNorms[pivot]);
448  std::swap(permutation[0], permutation[pivot]);
449  }
450  }
451 
452  qrHouseholderStepImpl(0, r, rhs, householder);
453 
454  MultiArrayIndex rank = 1;
455  NormType maxApproxSingularValue = norm(r(0,0)),
456  minApproxSingularValue = maxApproxSingularValue;
457 
458  double tolerance = (epsilon == 0.0)
459  ? m*maxApproxSingularValue*NumericTraits<T>::epsilon()
460  : epsilon;
461 
462  bool simpleSingularValueApproximation = (n < 4);
463  Matrix<T> zmax, zmin;
464  if(minApproxSingularValue <= tolerance)
465  {
466  rank = 0;
467  pivoting = false;
468  simpleSingularValueApproximation = true;
469  }
470  if(!simpleSingularValueApproximation)
471  {
472  zmax.reshape(Shape(m,1));
473  zmin.reshape(Shape(m,1));
474  zmax(0,0) = r(0,0);
475  zmin(0,0) = 1.0 / r(0,0);
476  }
477 
478  for(MultiArrayIndex k=1; k<maxRank; ++k)
479  {
480  if(pivoting)
481  {
482  for(MultiArrayIndex l=k; l<n; ++l)
483  columnSquaredNorms[l] -= squaredNorm(r(k, l));
484  int pivot = k + argMax(rowVector(columnSquaredNorms, Shape(0,k), n));
485  if(pivot != (int)k)
486  {
487  columnVector(r, k).swapData(columnVector(r, pivot));
488  std::swap(columnSquaredNorms[k], columnSquaredNorms[pivot]);
489  std::swap(permutation[k], permutation[pivot]);
490  }
491  }
492 
493  qrHouseholderStepImpl(k, r, rhs, householder);
494 
495  if(simpleSingularValueApproximation)
496  {
497  NormType nv = norm(r(k,k));
498  maxApproxSingularValue = std::max(nv, maxApproxSingularValue);
499  minApproxSingularValue = std::min(nv, minApproxSingularValue);
500  }
501  else
502  {
503  incrementalMaxSingularValueApproximation(columnVector(r, Shape(0,k),k+1), zmax, maxApproxSingularValue);
504  incrementalMinSingularValueApproximation(columnVector(r, Shape(0,k),k+1), zmin, minApproxSingularValue, tolerance);
505  }
506 
507 #if 0
508  Matrix<T> u(k+1,k+1), s(k+1, 1), v(k+1,k+1);
509  singularValueDecomposition(r.subarray(Shape(0,0), Shape(k+1,k+1)), u, s, v);
510  std::cerr << "estimate, svd " << k << ": " << minApproxSingularValue << " " << s(k,0) << "\n";
511 #endif
512 
513  if(epsilon == 0.0)
514  tolerance = m*maxApproxSingularValue*NumericTraits<T>::epsilon();
515 
516  if(minApproxSingularValue > tolerance)
517  ++rank;
518  else
519  pivoting = false; // matrix doesn't have full rank, triangulize the rest without pivoting
520  }
521  return (unsigned int)rank;
522 }
523 
524 template <class T, class C1, class C2>
525 unsigned int
526 qrTransformToUpperTriangular(MultiArrayView<2, T, C1> & r, MultiArrayView<2, T, C2> & rhs,
527  ArrayVector<MultiArrayIndex> & permutation, double epsilon = 0.0)
528 {
529  Matrix<T> dontStoreHouseholderVectors; // intentionally empty
530  return qrTransformToTriangularImpl(r, rhs, dontStoreHouseholderVectors, permutation, epsilon);
531 }
532 
533 // QR algorithm with optional row pivoting
534 template <class T, class C1, class C2, class C3>
535 unsigned int
536 qrTransformToLowerTriangular(MultiArrayView<2, T, C1> & r, MultiArrayView<2, T, C2> & rhs, MultiArrayView<2, T, C3> & householderMatrix,
537  double epsilon = 0.0)
538 {
539  ArrayVector<MultiArrayIndex> permutation((unsigned int)rowCount(rhs));
540  for(MultiArrayIndex k=0; k<(MultiArrayIndex)permutation.size(); ++k)
541  permutation[k] = k;
542  Matrix<T> dontTransformRHS; // intentionally empty
543  MultiArrayView<2, T, StridedArrayTag> rt = transpose(r),
544  ht = transpose(householderMatrix);
545  unsigned int rank = qrTransformToTriangularImpl(rt, dontTransformRHS, ht, permutation, epsilon);
546 
547  // apply row permutation to RHS
548  Matrix<T> tempRHS(rhs);
549  for(MultiArrayIndex k=0; k<(MultiArrayIndex)permutation.size(); ++k)
550  rowVector(rhs, k) = rowVector(tempRHS, permutation[k]);
551  return rank;
552 }
553 
554 // QR algorithm without column pivoting
555 template <class T, class C1, class C2>
556 inline bool
557 qrTransformToUpperTriangular(MultiArrayView<2, T, C1> & r, MultiArrayView<2, T, C2> & rhs,
558  double epsilon = 0.0)
559 {
560  ArrayVector<MultiArrayIndex> noPivoting; // intentionally empty
561 
562  return (qrTransformToUpperTriangular(r, rhs, noPivoting, epsilon) ==
563  (unsigned int)columnCount(r));
564 }
565 
566 // QR algorithm without row pivoting
567 template <class T, class C1, class C2>
568 inline bool
569 qrTransformToLowerTriangular(MultiArrayView<2, T, C1> & r, MultiArrayView<2, T, C2> & householder,
570  double epsilon = 0.0)
571 {
572  Matrix<T> noPivoting; // intentionally empty
573 
574  return (qrTransformToLowerTriangular(r, noPivoting, householder, epsilon) ==
575  (unsigned int)rowCount(r));
576 }
577 
578 // restore ordering of result vector elements after QR solution with column pivoting
579 template <class T, class C1, class C2, class Permutation>
580 void inverseRowPermutation(MultiArrayView<2, T, C1> &permuted, MultiArrayView<2, T, C2> &res,
581  Permutation const & permutation)
582 {
583  for(MultiArrayIndex k=0; k<columnCount(permuted); ++k)
584  for(MultiArrayIndex l=0; l<rowCount(permuted); ++l)
585  res(permutation[l], k) = permuted(l,k);
586 }
587 
588 template <class T, class C1, class C2>
589 void applyHouseholderColumnReflections(MultiArrayView<2, T, C1> const &householder, MultiArrayView<2, T, C2> &res)
590 {
591  typedef typename Matrix<T>::difference_type Shape;
592  MultiArrayIndex n = rowCount(householder);
593  MultiArrayIndex m = columnCount(householder);
594  MultiArrayIndex rhsCount = columnCount(res);
595 
596  for(int k = m-1; k >= 0; --k)
597  {
598  MultiArrayView<2, T, C1> u = columnVector(householder, Shape(k,k), n);
599  for(MultiArrayIndex l=0; l<rhsCount; ++l)
600  columnVector(res, Shape(k,l), n) -= dot(columnVector(res, Shape(k,l), n), u) * u;
601  }
602 }
603 
604 } // namespace detail
605 
606 template <class T, class C1, class C2, class C3>
607 unsigned int
608 linearSolveQRReplace(MultiArrayView<2, T, C1> &A, MultiArrayView<2, T, C2> &b,
609  MultiArrayView<2, T, C3> & res,
610  double epsilon = 0.0)
611 {
612  typedef typename Matrix<T>::difference_type Shape;
613 
615  MultiArrayIndex m = rowCount(A);
616  MultiArrayIndex rhsCount = columnCount(res);
617  MultiArrayIndex rank = std::min(m,n);
618  Shape ul(MultiArrayIndex(0), MultiArrayIndex(0));
619 
620 
621  vigra_precondition(rhsCount == columnCount(b),
622  "linearSolveQR(): RHS and solution must have the same number of columns.");
623  vigra_precondition(m == rowCount(b),
624  "linearSolveQR(): Coefficient matrix and RHS must have the same number of rows.");
625  vigra_precondition(n == rowCount(res),
626  "linearSolveQR(): Mismatch between column count of coefficient matrix and row count of solution.");
627  vigra_precondition(epsilon >= 0.0,
628  "linearSolveQR(): 'epsilon' must be non-negative.");
629 
630  if(m < n)
631  {
632  // minimum norm solution of underdetermined system
633  Matrix<T> householderMatrix(n, m);
634  MultiArrayView<2, T, StridedArrayTag> ht = transpose(householderMatrix);
635  rank = (MultiArrayIndex)detail::qrTransformToLowerTriangular(A, b, ht, epsilon);
636  res.subarray(Shape(rank,0), Shape(n, rhsCount)).init(NumericTraits<T>::zero());
637  if(rank < m)
638  {
639  // system is also rank-deficient => compute minimum norm least squares solution
640  MultiArrayView<2, T, C1> Asub = A.subarray(ul, Shape(m,rank));
641  detail::qrTransformToUpperTriangular(Asub, b, epsilon);
642  linearSolveUpperTriangular(A.subarray(ul, Shape(rank,rank)),
643  b.subarray(ul, Shape(rank,rhsCount)),
644  res.subarray(ul, Shape(rank, rhsCount)));
645  }
646  else
647  {
648  // system has full rank => compute minimum norm solution
649  linearSolveLowerTriangular(A.subarray(ul, Shape(rank,rank)),
650  b.subarray(ul, Shape(rank, rhsCount)),
651  res.subarray(ul, Shape(rank, rhsCount)));
652  }
653  detail::applyHouseholderColumnReflections(householderMatrix.subarray(ul, Shape(n, rank)), res);
654  }
655  else
656  {
657  // solution of well-determined or overdetermined system
658  ArrayVector<MultiArrayIndex> permutation((unsigned int)n);
659  for(MultiArrayIndex k=0; k<n; ++k)
660  permutation[k] = k;
661 
662  rank = (MultiArrayIndex)detail::qrTransformToUpperTriangular(A, b, permutation, epsilon);
663 
664  Matrix<T> permutedSolution(n, rhsCount);
665  if(rank < n)
666  {
667  // system is rank-deficient => compute minimum norm solution
668  Matrix<T> householderMatrix(n, rank);
669  MultiArrayView<2, T, StridedArrayTag> ht = transpose(householderMatrix);
670  MultiArrayView<2, T, C1> Asub = A.subarray(ul, Shape(rank,n));
671  detail::qrTransformToLowerTriangular(Asub, ht, epsilon);
672  linearSolveLowerTriangular(A.subarray(ul, Shape(rank,rank)),
673  b.subarray(ul, Shape(rank, rhsCount)),
674  permutedSolution.subarray(ul, Shape(rank, rhsCount)));
675  detail::applyHouseholderColumnReflections(householderMatrix, permutedSolution);
676  }
677  else
678  {
679  // system has full rank => compute exact or least squares solution
680  linearSolveUpperTriangular(A.subarray(ul, Shape(rank,rank)),
681  b.subarray(ul, Shape(rank,rhsCount)),
682  permutedSolution);
683  }
684  detail::inverseRowPermutation(permutedSolution, res, permutation);
685  }
686  return (unsigned int)rank;
687 }
688 
689 template <class T, class C1, class C2, class C3>
690 unsigned int linearSolveQR(MultiArrayView<2, T, C1> const & A, MultiArrayView<2, T, C2> const & b,
691  MultiArrayView<2, T, C3> & res)
692 {
693  Matrix<T> r(A), rhs(b);
694  return linearSolveQRReplace(r, rhs, res);
695 }
696 
697 /** \defgroup MatrixAlgebra Advanced Matrix Algebra
698 
699  \brief Solution of linear systems, eigen systems, linear least squares etc.
700 
701  \ingroup LinearAlgebraModule
702  */
703 //@{
704  /** Create the inverse or pseudo-inverse of matrix \a v.
705 
706  If the matrix \a v is square, \a res must have the same shape and will contain the
707  inverse of \a v. If \a v is rectangular, \a res must have the transposed shape
708  of \a v. The inverse is then computed in the least-squares
709  sense, i.e. \a res will be the pseudo-inverse (Moore-Penrose inverse).
710  The function returns <tt>true</tt> upon success, and <tt>false</tt> if \a v
711  is not invertible (has not full rank). The inverse is computed by means of QR
712  decomposition. This function can be applied in-place.
713 
714  <b>\#include</b> <vigra/linear_solve.hxx> or<br>
715  <b>\#include</b> <vigra/linear_algebra.hxx><br>
716  Namespaces: vigra and vigra::linalg
717  */
718 template <class T, class C1, class C2>
720 {
721  typedef typename MultiArrayShape<2>::type Shape;
722 
723  const MultiArrayIndex n = columnCount(v);
724  const MultiArrayIndex m = rowCount(v);
725  vigra_precondition(n == rowCount(res) && m == columnCount(res),
726  "inverse(): shape of output matrix must be the transpose of the input matrix' shape.");
727 
728  if(m < n)
729  {
731  Matrix<T> r(vt.shape()), q(n, n);
732  if(!qrDecomposition(vt, q, r))
733  return false; // a didn't have full rank
734  linearSolveUpperTriangular(r.subarray(Shape(0,0), Shape(m,m)),
735  transpose(q).subarray(Shape(0,0), Shape(m,n)),
736  transpose(res));
737  }
738  else
739  {
740  Matrix<T> r(v.shape()), q(m, m);
741  if(!qrDecomposition(v, q, r))
742  return false; // a didn't have full rank
743  linearSolveUpperTriangular(r.subarray(Shape(0,0), Shape(n,n)),
744  transpose(q).subarray(Shape(0,0), Shape(n,m)),
745  res);
746  }
747  return true;
748 }
749 
750  /** Create the inverse or pseudo-inverse of matrix \a v.
751 
752  The result is returned as a temporary matrix. If the matrix \a v is square,
753  the result will have the same shape and contains the inverse of \a v.
754  If \a v is rectangular, the result will have the transposed shape of \a v.
755  The inverse is then computed in the least-squares
756  sense, i.e. \a res will be the pseudo-inverse (Moore-Penrose inverse).
757  The inverse is computed by means of QR decomposition. If \a v
758  is not invertible, <tt>vigra::PreconditionViolation</tt> exception is thrown.
759  Usage:
760 
761  \code
762  vigra::Matrix<double> v(n, n);
763  v = ...;
764 
765  vigra::Matrix<double> m = inverse(v);
766  \endcode
767 
768  <b>\#include</b> <vigra/linear_solve.hxx> or<br>
769  <b>\#include</b> <vigra/linear_algebra.hxx><br>
770  Namespaces: vigra and vigra::linalg
771  */
772 template <class T, class C>
773 TemporaryMatrix<T> inverse(const MultiArrayView<2, T, C> &v)
774 {
775  TemporaryMatrix<T> ret(columnCount(v), rowCount(v)); // transpose shape
776  vigra_precondition(inverse(v, ret),
777  "inverse(): matrix is not invertible.");
778  return ret;
779 }
780 
781 template <class T>
782 TemporaryMatrix<T> inverse(const TemporaryMatrix<T> &v)
783 {
784  if(columnCount(v) == rowCount(v))
785  {
786  vigra_precondition(inverse(v, const_cast<TemporaryMatrix<T> &>(v)),
787  "inverse(): matrix is not invertible.");
788  return v;
789  }
790  else
791  {
792  TemporaryMatrix<T> ret(columnCount(v), rowCount(v)); // transpose shape
793  vigra_precondition(inverse(v, ret),
794  "inverse(): matrix is not invertible.");
795  return ret;
796  }
797 }
798 
799  /** Compute the determinant of a square matrix.
800 
801  \a method must be one of the following:
802  <DL>
803  <DT>"Cholesky"<DD> Compute the solution by means of Cholesky decomposition. This
804  method is faster than "LU", but requires the matrix \a a
805  to be symmetric positive definite. If this is
806  not the case, a <tt>ContractViolation</tt> exception is thrown.
807 
808  <DT>"LU"<DD> (default) Compute the solution by means of LU decomposition.
809  </DL>
810 
811  <b>\#include</b> <vigra/linear_solve.hxx> or<br>
812  <b>\#include</b> <vigra/linear_algebra.hxx><br>
813  Namespaces: vigra and vigra::linalg
814  */
815 template <class T, class C1>
816 T determinant(MultiArrayView<2, T, C1> const & a, std::string method = "LU")
817 {
819  vigra_precondition(rowCount(a) == n,
820  "determinant(): Square matrix required.");
821 
822  for(unsigned int k=0; k<method.size(); ++k)
823  method[k] = tolower(method[k]);
824 
825  if(n == 1)
826  return a(0,0);
827  if(n == 2)
828  return a(0,0)*a(1,1) - a(0,1)*a(1,0);
829  if(method == "lu")
830  {
831  return detail::determinantByLUDecomposition(a);
832  }
833  else if(method == "cholesky")
834  {
835  Matrix<T> L(a.shape());
836  vigra_precondition(choleskyDecomposition(a, L),
837  "determinant(): Cholesky method requires symmetric positive definite matrix.");
838  T det = L(0,0);
839  for(MultiArrayIndex k=1; k<n; ++k)
840  det *= L(k,k);
841  return sq(det);
842  }
843  else
844  {
845  vigra_precondition(false, "determinant(): Unknown solution method.");
846  }
847  return T();
848 }
849 
850  /** Compute the logarithm of the determinant of a symmetric positive definite matrix.
851 
852  This is useful to avoid multiplication of very large numbers in big matrices.
853  It is implemented by means of Cholesky decomposition.
854 
855  <b>\#include</b> <vigra/linear_solve.hxx> or<br>
856  <b>\#include</b> <vigra/linear_algebra.hxx><br>
857  Namespaces: vigra and vigra::linalg
858  */
859 template <class T, class C1>
861 {
863  vigra_precondition(rowCount(a) == n,
864  "logDeterminant(): Square matrix required.");
865  if(n == 1)
866  {
867  vigra_precondition(a(0,0) > 0.0,
868  "logDeterminant(): Matrix not positive definite.");
869  return std::log(a(0,0));
870  }
871  if(n == 2)
872  {
873  T det = a(0,0)*a(1,1) - a(0,1)*a(1,0);
874  vigra_precondition(det > 0.0,
875  "logDeterminant(): Matrix not positive definite.");
876  return std::log(det);
877  }
878  else
879  {
880  Matrix<T> L(a.shape());
881  vigra_precondition(choleskyDecomposition(a, L),
882  "logDeterminant(): Matrix not positive definite.");
883  T logdet = std::log(L(0,0));
884  for(MultiArrayIndex k=1; k<n; ++k)
885  logdet += std::log(L(k,k)); // L(k,k) is guaranteed to be positive
886  return 2.0*logdet;
887  }
888 }
889 
890  /** Cholesky decomposition.
891 
892  \a A must be a symmetric positive definite matrix, and \a L will be a lower
893  triangular matrix, such that (up to round-off errors):
894 
895  \code
896  A == L * transpose(L);
897  \endcode
898 
899  This implementation cannot be applied in-place, i.e. <tt>&L == &A</tt> is an error.
900  If \a A is not symmetric, a <tt>ContractViolation</tt> exception is thrown. If it
901  is not positive definite, the function returns <tt>false</tt>.
902 
903  <b>\#include</b> <vigra/linear_solve.hxx> or<br>
904  <b>\#include</b> <vigra/linear_algebra.hxx><br>
905  Namespaces: vigra and vigra::linalg
906  */
907 template <class T, class C1, class C2>
910 {
912  vigra_precondition(rowCount(A) == n,
913  "choleskyDecomposition(): Input matrix must be square.");
914  vigra_precondition(n == columnCount(L) && n == rowCount(L),
915  "choleskyDecomposition(): Output matrix must have same shape as input matrix.");
916  vigra_precondition(isSymmetric(A),
917  "choleskyDecomposition(): Input matrix must be symmetric.");
918 
919  for (MultiArrayIndex j = 0; j < n; ++j)
920  {
921  T d(0.0);
922  for (MultiArrayIndex k = 0; k < j; ++k)
923  {
924  T s(0.0);
925  for (MultiArrayIndex i = 0; i < k; ++i)
926  {
927  s += L(k, i)*L(j, i);
928  }
929  L(j, k) = s = (A(j, k) - s)/L(k, k);
930  d = d + s*s;
931  }
932  d = A(j, j) - d;
933  if(d <= 0.0)
934  return false; // A is not positive definite
935  L(j, j) = std::sqrt(d);
936  for (MultiArrayIndex k = j+1; k < n; ++k)
937  {
938  L(j, k) = 0.0;
939  }
940  }
941  return true;
942 }
943 
944  /** QR decomposition.
945 
946  \a a contains the original matrix, results are returned in \a q and \a r, where
947  \a q is a orthogonal matrix, and \a r is an upper triangular matrix, such that
948  (up to round-off errors):
949 
950  \code
951  a == q * r;
952  \endcode
953 
954  If \a a doesn't have full rank, the function returns <tt>false</tt>.
955  The decomposition is computed by householder transformations. It can be applied in-place,
956  i.e. <tt>&a == &q</tt> or <tt>&a == &r</tt> are allowed.
957 
958  <b>\#include</b> <vigra/linear_solve.hxx> or<br>
959  <b>\#include</b> <vigra/linear_algebra.hxx><br>
960  Namespaces: vigra and vigra::linalg
961  */
962 template <class T, class C1, class C2, class C3>
965  double epsilon = 0.0)
966 {
967  const MultiArrayIndex m = rowCount(a);
968  const MultiArrayIndex n = columnCount(a);
969  vigra_precondition(n == columnCount(r) && m == rowCount(r) &&
970  m == columnCount(q) && m == rowCount(q),
971  "qrDecomposition(): Matrix shape mismatch.");
972 
973  q = identityMatrix<T>(m);
975  r = a;
976  ArrayVector<MultiArrayIndex> noPivoting; // intentionally empty
977  return ((MultiArrayIndex)detail::qrTransformToUpperTriangular(r, tq, noPivoting, epsilon) == std::min(m,n));
978 }
979 
980  /** Deprecated, use \ref linearSolveUpperTriangular().
981  */
982 template <class T, class C1, class C2, class C3>
983 inline
986 {
987  return linearSolveUpperTriangular(r, b, x);
988 }
989 
990  /** Solve a linear system with upper-triangular coefficient matrix.
991 
992  The square matrix \a r must be an upper-triangular coefficient matrix as can,
993  for example, be obtained by means of QR decomposition. If \a r doesn't have full rank
994  the function fails and returns <tt>false</tt>, otherwise it returns <tt>true</tt>. The
995  lower triangular part of matrix \a r will not be touched, so it doesn't need to contain zeros.
996 
997  The column vectors of matrix \a b are the right-hand sides of the equation (several equations
998  with the same coefficients can thus be solved in one go). The result is returned
999  int \a x, whose columns contain the solutions for the corresponding
1000  columns of \a b. This implementation can be applied in-place, i.e. <tt>&b == &x</tt> is allowed.
1001  The following size requirements apply:
1002 
1003  \code
1004  rowCount(r) == columnCount(r);
1005  rowCount(r) == rowCount(b);
1006  columnCount(r) == rowCount(x);
1007  columnCount(b) == columnCount(x);
1008  \endcode
1009 
1010  <b>\#include</b> <vigra/linear_solve.hxx> or<br>
1011  <b>\#include</b> <vigra/linear_algebra.hxx><br>
1012  Namespaces: vigra and vigra::linalg
1013  */
1014 template <class T, class C1, class C2, class C3>
1017 {
1018  typedef MultiArrayShape<2>::type Shape;
1019  MultiArrayIndex m = rowCount(r);
1020  MultiArrayIndex rhsCount = columnCount(b);
1021  vigra_precondition(m == columnCount(r),
1022  "linearSolveUpperTriangular(): square coefficient matrix required.");
1023  vigra_precondition(m == rowCount(b) && m == rowCount(x) && rhsCount == columnCount(x),
1024  "linearSolveUpperTriangular(): matrix shape mismatch.");
1025 
1026  for(MultiArrayIndex k = 0; k < rhsCount; ++k)
1027  {
1028  for(int i=m-1; i>=0; --i)
1029  {
1030  if(r(i,i) == NumericTraits<T>::zero())
1031  return false; // r doesn't have full rank
1032  T sum = b(i, k);
1033  for(MultiArrayIndex j=i+1; j<m; ++j)
1034  sum -= r(i, j) * x(j, k);
1035  x(i, k) = sum / r(i, i);
1036  }
1037  }
1038  return true;
1039 }
1040 
1041  /** Solve a linear system with lower-triangular coefficient matrix.
1042 
1043  The square matrix \a l must be a lower-triangular coefficient matrix. If \a l
1044  doesn't have full rank the function fails and returns <tt>false</tt>,
1045  otherwise it returns <tt>true</tt>. The upper triangular part of matrix \a l will not be touched,
1046  so it doesn't need to contain zeros.
1047 
1048  The column vectors of matrix \a b are the right-hand sides of the equation (several equations
1049  with the same coefficients can thus be solved in one go). The result is returned
1050  in \a x, whose columns contain the solutions for the corresponding
1051  columns of \a b. This implementation can be applied in-place, i.e. <tt>&b == &x</tt> is allowed.
1052  The following size requirements apply:
1053 
1054  \code
1055  rowCount(l) == columnCount(l);
1056  rowCount(l) == rowCount(b);
1057  columnCount(l) == rowCount(x);
1058  columnCount(b) == columnCount(x);
1059  \endcode
1060 
1061  <b>\#include</b> <vigra/linear_solve.hxx> or<br>
1062  <b>\#include</b> <vigra/linear_algebra.hxx><br>
1063  Namespaces: vigra and vigra::linalg
1064  */
1065 template <class T, class C1, class C2, class C3>
1068 {
1071  vigra_precondition(m == rowCount(l),
1072  "linearSolveLowerTriangular(): square coefficient matrix required.");
1073  vigra_precondition(m == rowCount(b) && m == rowCount(x) && n == columnCount(x),
1074  "linearSolveLowerTriangular(): matrix shape mismatch.");
1075 
1076  for(MultiArrayIndex k = 0; k < n; ++k)
1077  {
1078  for(MultiArrayIndex i=0; i<m; ++i)
1079  {
1080  if(l(i,i) == NumericTraits<T>::zero())
1081  return false; // l doesn't have full rank
1082  T sum = b(i, k);
1083  for(MultiArrayIndex j=0; j<i; ++j)
1084  sum -= l(i, j) * x(j, k);
1085  x(i, k) = sum / l(i, i);
1086  }
1087  }
1088  return true;
1089 }
1090 
1091 
1092  /** Solve a linear system when the Cholesky decomposition of the left hand side is given.
1093 
1094  The square matrix \a L must be a lower-triangular matrix resulting from Cholesky
1095  decomposition of some positive definite coefficient matrix.
1096 
1097  The column vectors of matrix \a b are the right-hand sides of the equation (several equations
1098  with the same matrix \a L can thus be solved in one go). The result is returned
1099  in \a x, whose columns contain the solutions for the corresponding
1100  columns of \a b. This implementation can be applied in-place, i.e. <tt>&b == &x</tt> is allowed.
1101  The following size requirements apply:
1102 
1103  \code
1104  rowCount(L) == columnCount(L);
1105  rowCount(L) == rowCount(b);
1106  columnCount(L) == rowCount(x);
1107  columnCount(b) == columnCount(x);
1108  \endcode
1109 
1110  <b>\#include</b> <vigra/linear_solve.hxx> or<br>
1111  <b>\#include</b> <vigra/linear_algebra.hxx><br>
1112  Namespaces: vigra and vigra::linalg
1113  */
1114 template <class T, class C1, class C2, class C3>
1115 inline
1117 {
1118  /* Solve L * y = b */
1119  linearSolveLowerTriangular(L, b, x);
1120  /* Solve L^T * x = y */
1122 }
1123 
1124  /** Solve a linear system.
1125 
1126  \a A is the coefficient matrix, and the column vectors
1127  in \a b are the right-hand sides of the equation (so, several equations
1128  with the same coefficients can be solved in one go). The result is returned
1129  in \a res, whose columns contain the solutions for the corresponding
1130  columns of \a b. The number of columns of \a A must equal the number of rows of
1131  both \a b and \a res, and the number of columns of \a b and \a res must match.
1132 
1133  \a method must be one of the following:
1134  <DL>
1135  <DT>"Cholesky"<DD> Compute the solution by means of Cholesky decomposition. The
1136  coefficient matrix \a A must by symmetric positive definite. If
1137  this is not the case, the function returns <tt>false</tt>.
1138 
1139  <DT>"QR"<DD> (default) Compute the solution by means of QR decomposition. The
1140  coefficient matrix \a A can be square or rectangular. In the latter case,
1141  it must have more rows than columns, and the solution will be computed in the
1142  least squares sense. If \a A doesn't have full rank, the function
1143  returns <tt>false</tt>.
1144 
1145  <DT>"SVD"<DD> Compute the solution by means of singular value decomposition. The
1146  coefficient matrix \a A can be square or rectangular. In the latter case,
1147  it must have more rows than columns, and the solution will be computed in the
1148  least squares sense. If \a A doesn't have full rank, the function
1149  returns <tt>false</tt>.
1150 
1151  <DT>"NE"<DD> Compute the solution by means of the normal equations, i.e. by applying Cholesky
1152  decomposition to the equivalent problem <tt>A'*A*x = A'*b</tt>. This only makes sense
1153  when the equation is to be solved in the least squares sense, i.e. when \a A is a
1154  rectangular matrix with more rows than columns. If \a A doesn't have full column rank,
1155  the function returns <tt>false</tt>.
1156  </DL>
1157 
1158  This function can be applied in-place, i.e. <tt>&b == &res</tt> or <tt>&A == &res</tt> are allowed
1159  (provided they have the required shapes).
1160 
1161  The following size requirements apply:
1162 
1163  \code
1164  rowCount(r) == rowCount(b);
1165  columnCount(r) == rowCount(x);
1166  columnCount(b) == columnCount(x);
1167  \endcode
1168 
1169  <b>\#include</b> <vigra/linear_solve.hxx> or<br>
1170  <b>\#include</b> <vigra/linear_algebra.hxx><br>
1171  Namespaces: vigra and vigra::linalg
1172  */
1173 template <class T, class C1, class C2, class C3>
1175  MultiArrayView<2, T, C3> & res, std::string method = "QR")
1176 {
1177  typedef typename Matrix<T>::difference_type Shape;
1178  typedef typename Matrix<T>::view_type SubMatrix;
1179 
1180  const MultiArrayIndex n = columnCount(A);
1181  const MultiArrayIndex m = rowCount(A);
1182 
1183  vigra_precondition(n <= m,
1184  "linearSolve(): Coefficient matrix A must have at least as many rows as columns.");
1185  vigra_precondition(n == rowCount(res) &&
1186  m == rowCount(b) && columnCount(b) == columnCount(res),
1187  "linearSolve(): matrix shape mismatch.");
1188 
1189  for(unsigned int k=0; k<method.size(); ++k)
1190  method[k] = (std::string::value_type)tolower(method[k]);
1191 
1192  if(method == "cholesky")
1193  {
1194  vigra_precondition(columnCount(A) == rowCount(A),
1195  "linearSolve(): Cholesky method requires square coefficient matrix.");
1196  Matrix<T> L(A.shape());
1197  if(!choleskyDecomposition(A, L))
1198  return false; // false if A wasn't symmetric positive definite
1199  choleskySolve(L, b, res);
1200  }
1201  else if(method == "qr")
1202  {
1203  return (MultiArrayIndex)linearSolveQR(A, b, res) == n;
1204  }
1205  else if(method == "ne")
1206  {
1207  return linearSolve(transpose(A)*A, transpose(A)*b, res, "Cholesky");
1208  }
1209  else if(method == "svd")
1210  {
1211  MultiArrayIndex rhsCount = columnCount(b);
1212  Matrix<T> u(A.shape()), s(n, 1), v(n, n);
1213 
1215 
1216  Matrix<T> t = transpose(u)*b;
1217  for(MultiArrayIndex l=0; l<rhsCount; ++l)
1218  {
1219  for(MultiArrayIndex k=0; k<rank; ++k)
1220  t(k,l) /= s(k,0);
1221  for(MultiArrayIndex k=rank; k<n; ++k)
1222  t(k,l) = NumericTraits<T>::zero();
1223  }
1224  res = v*t;
1225 
1226  return (rank == n);
1227  }
1228  else
1229  {
1230  vigra_precondition(false, "linearSolve(): Unknown solution method.");
1231  }
1232  return true;
1233 }
1234 
1235 //@}
1236 
1237 } // namespace linalg
1238 
1239 using linalg::inverse;
1240 using linalg::determinant;
1242 using linalg::linearSolve;
1243 using linalg::choleskySolve;
1248 
1249 } // namespace vigra
1250 
1251 
1252 #endif // VIGRA_LINEAR_SOLVE_HXX

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
vigra 1.8.0 (Wed Sep 26 2012)