[ 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, it must have more rows than columns, and \a res
708  must have the transposed shape 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> <<a href="linear__solve_8hxx-source.html">vigra/linear_solve.hxx</a>> or<br>
715  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
716  Namespaces: vigra and vigra::linalg
717  */
718 template <class T, class C1, class C2>
720 {
721  const MultiArrayIndex n = columnCount(v);
722  vigra_precondition(n <= rowCount(v),
723  "inverse(): input matrix must have at least as many rows as columns.");
724  vigra_precondition(n == rowCount(res) && rowCount(v) == columnCount(res),
725  "inverse(): shape of output matrix must be the transpose of the input matrix' shape.");
726 
727  Matrix<T> r(v.shape()), q(n, n);
728  if(!qrDecomposition(v, q, r))
729  return false; // a didn't have full rank
731  return true;
732 }
733 
734  /** Create the inverse or pseudo-inverse of matrix \a v.
735 
736  The result is returned as a temporary matrix. If the matrix \a v is square,
737  the result will have the same shape and contains the inverse of \a v.
738  If \a v is rectangular, it must have more rows than columns, and the result will
739  have the transposed shape of \a v. The inverse is then computed in the least-squares
740  sense, i.e. \a res will be the pseudo-inverse (Moore-Penrose inverse).
741  The inverse is computed by means of QR decomposition. If \a v
742  is not invertible, <tt>vigra::PreconditionViolation</tt> exception is thrown.
743  Usage:
744 
745  \code
746  vigra::Matrix<double> v(n, n);
747  v = ...;
748 
749  vigra::Matrix<double> m = inverse(v);
750  \endcode
751 
752  <b>\#include</b> <<a href="linear__solve_8hxx-source.html">vigra/linear_solve.hxx</a>> or<br>
753  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
754  Namespaces: vigra and vigra::linalg
755  */
756 template <class T, class C>
757 TemporaryMatrix<T> inverse(const MultiArrayView<2, T, C> &v)
758 {
759  TemporaryMatrix<T> ret(columnCount(v), rowCount(v)); // transpose shape
760  vigra_precondition(inverse(v, ret),
761  "inverse(): matrix is not invertible.");
762  return ret;
763 }
764 
765 template <class T>
766 TemporaryMatrix<T> inverse(const TemporaryMatrix<T> &v)
767 {
768  if(columnCount(v) == rowCount(v))
769  {
770  vigra_precondition(inverse(v, const_cast<TemporaryMatrix<T> &>(v)),
771  "inverse(): matrix is not invertible.");
772  return v;
773  }
774  else
775  {
776  TemporaryMatrix<T> ret(columnCount(v), rowCount(v)); // transpose shape
777  vigra_precondition(inverse(v, ret),
778  "inverse(): matrix is not invertible.");
779  return ret;
780  }
781 }
782 
783  /** Compute the determinant of a square matrix.
784 
785  \a method must be one of the following:
786  <DL>
787  <DT>"Cholesky"<DD> Compute the solution by means of Cholesky decomposition. This
788  method is faster than "LU", but requires the matrix \a a
789  to be symmetric positive definite. If this is
790  not the case, a <tt>ContractViolation</tt> exception is thrown.
791 
792  <DT>"LU"<DD> (default) Compute the solution by means of LU decomposition.
793  </DL>
794 
795  <b>\#include</b> <<a href="linear__solve_8hxx-source.html">vigra/linear_solve.hxx</a>> or<br>
796  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
797  Namespaces: vigra and vigra::linalg
798  */
799 template <class T, class C1>
800 T determinant(MultiArrayView<2, T, C1> const & a, std::string method = "LU")
801 {
803  vigra_precondition(rowCount(a) == n,
804  "determinant(): Square matrix required.");
805 
806  for(unsigned int k=0; k<method.size(); ++k)
807  method[k] = tolower(method[k]);
808 
809  if(n == 1)
810  return a(0,0);
811  if(n == 2)
812  return a(0,0)*a(1,1) - a(0,1)*a(1,0);
813  if(method == "lu")
814  {
815  return detail::determinantByLUDecomposition(a);
816  }
817  else if(method == "cholesky")
818  {
819  Matrix<T> L(a.shape());
820  vigra_precondition(choleskyDecomposition(a, L),
821  "determinant(): Cholesky method requires symmetric positive definite matrix.");
822  T det = L(0,0);
823  for(MultiArrayIndex k=1; k<n; ++k)
824  det *= L(k,k);
825  return sq(det);
826  }
827  else
828  {
829  vigra_precondition(false, "determinant(): Unknown solution method.");
830  }
831  return T();
832 }
833 
834  /** Compute the logarithm of the determinant of a symmetric positive definite matrix.
835 
836  This is useful to avoid multiplication of very large numbers in big matrices.
837  It is implemented by means of Cholesky decomposition.
838 
839  <b>\#include</b> <<a href="linear__solve_8hxx-source.html">vigra/linear_solve.hxx</a>> or<br>
840  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
841  Namespaces: vigra and vigra::linalg
842  */
843 template <class T, class C1>
845 {
847  vigra_precondition(rowCount(a) == n,
848  "logDeterminant(): Square matrix required.");
849  if(n == 1)
850  {
851  vigra_precondition(a(0,0) > 0.0,
852  "logDeterminant(): Matrix not positive definite.");
853  return std::log(a(0,0));
854  }
855  if(n == 2)
856  {
857  T det = a(0,0)*a(1,1) - a(0,1)*a(1,0);
858  vigra_precondition(det > 0.0,
859  "logDeterminant(): Matrix not positive definite.");
860  return std::log(det);
861  }
862  else
863  {
864  Matrix<T> L(a.shape());
865  vigra_precondition(choleskyDecomposition(a, L),
866  "logDeterminant(): Matrix not positive definite.");
867  T logdet = std::log(L(0,0));
868  for(MultiArrayIndex k=1; k<n; ++k)
869  logdet += std::log(L(k,k)); // L(k,k) is guaranteed to be positive
870  return 2.0*logdet;
871  }
872 }
873 
874  /** Cholesky decomposition.
875 
876  \a A must be a symmetric positive definite matrix, and \a L will be a lower
877  triangular matrix, such that (up to round-off errors):
878 
879  \code
880  A == L * transpose(L);
881  \endcode
882 
883  This implementation cannot be applied in-place, i.e. <tt>&L == &A</tt> is an error.
884  If \a A is not symmetric, a <tt>ContractViolation</tt> exception is thrown. If it
885  is not positive definite, the function returns <tt>false</tt>.
886 
887  <b>\#include</b> <<a href="linear__solve_8hxx-source.html">vigra/linear_solve.hxx</a>> or<br>
888  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
889  Namespaces: vigra and vigra::linalg
890  */
891 template <class T, class C1, class C2>
894 {
896  vigra_precondition(rowCount(A) == n,
897  "choleskyDecomposition(): Input matrix must be square.");
898  vigra_precondition(n == columnCount(L) && n == rowCount(L),
899  "choleskyDecomposition(): Output matrix must have same shape as input matrix.");
900  vigra_precondition(isSymmetric(A),
901  "choleskyDecomposition(): Input matrix must be symmetric.");
902 
903  for (MultiArrayIndex j = 0; j < n; ++j)
904  {
905  T d(0.0);
906  for (MultiArrayIndex k = 0; k < j; ++k)
907  {
908  T s(0.0);
909  for (MultiArrayIndex i = 0; i < k; ++i)
910  {
911  s += L(k, i)*L(j, i);
912  }
913  L(j, k) = s = (A(j, k) - s)/L(k, k);
914  d = d + s*s;
915  }
916  d = A(j, j) - d;
917  if(d <= 0.0)
918  return false; // A is not positive definite
919  L(j, j) = std::sqrt(d);
920  for (MultiArrayIndex k = j+1; k < n; ++k)
921  {
922  L(j, k) = 0.0;
923  }
924  }
925  return true;
926 }
927 
928  /** QR decomposition.
929 
930  \a a contains the original matrix, results are returned in \a q and \a r, where
931  \a q is a orthogonal matrix, and \a r is an upper triangular matrix, such that
932  (up to round-off errors):
933 
934  \code
935  a == q * r;
936  \endcode
937 
938  If \a a dosn't have full rank, the function returns <tt>false</tt>.
939  The decomposition is computed by householder transformations. It can be applied in-place,
940  i.e. <tt>&a == &q</tt> or <tt>&a == &r</tt> are allowed.
941 
942  <b>\#include</b> <<a href="linear__solve_8hxx-source.html">vigra/linear_solve.hxx</a>> or<br>
943  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
944  Namespaces: vigra and vigra::linalg
945  */
946 template <class T, class C1, class C2, class C3>
949  double epsilon = 0.0)
950 {
951  const MultiArrayIndex m = rowCount(a);
952  const MultiArrayIndex n = columnCount(a);
953  vigra_precondition(n == columnCount(r) && m == rowCount(r) &&
954  m == columnCount(q) && m == rowCount(q),
955  "qrDecomposition(): Matrix shape mismatch.");
956 
957  q = identityMatrix<T>(m);
959  r = a;
960  ArrayVector<MultiArrayIndex> noPivoting; // intentionally empty
961  return ((MultiArrayIndex)detail::qrTransformToUpperTriangular(r, tq, noPivoting, epsilon) == std::min(m,n));
962 }
963 
964  /** Deprecated, use \ref linearSolveUpperTriangular().
965  */
966 template <class T, class C1, class C2, class C3>
967 inline
970 {
971  return linearSolveUpperTriangular(r, b, x);
972 }
973 
974  /** Solve a linear system with upper-triangular coefficient matrix.
975 
976  The square matrix \a r must be an upper-triangular coefficient matrix as can,
977  for example, be obtained by means of QR decomposition. If \a r doesn't have full rank
978  the function fails and returns <tt>false</tt>, otherwise it returns <tt>true</tt>. The
979  lower triangular part of matrix \a r will not be touched, so it doesn't need to contain zeros.
980 
981  The column vectors of matrix \a b are the right-hand sides of the equation (several equations
982  with the same coefficients can thus be solved in one go). The result is returned
983  int \a x, whose columns contain the solutions for the corresponding
984  columns of \a b. This implementation can be applied in-place, i.e. <tt>&b == &x</tt> is allowed.
985  The following size requirements apply:
986 
987  \code
988  rowCount(r) == columnCount(r);
989  rowCount(r) == rowCount(b);
990  columnCount(r) == rowCount(x);
991  columnCount(b) == columnCount(x);
992  \endcode
993 
994  <b>\#include</b> <<a href="linear__solve_8hxx-source.html">vigra/linear_solve.hxx</a>> or<br>
995  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
996  Namespaces: vigra and vigra::linalg
997  */
998 template <class T, class C1, class C2, class C3>
1001 {
1002  typedef MultiArrayShape<2>::type Shape;
1003  MultiArrayIndex m = rowCount(r);
1004  MultiArrayIndex rhsCount = columnCount(b);
1005  vigra_precondition(m == columnCount(r),
1006  "linearSolveUpperTriangular(): square coefficient matrix required.");
1007  vigra_precondition(m == rowCount(b) && m == rowCount(x) && rhsCount == columnCount(x),
1008  "linearSolveUpperTriangular(): matrix shape mismatch.");
1009 
1010  for(MultiArrayIndex k = 0; k < rhsCount; ++k)
1011  {
1012  for(int i=m-1; i>=0; --i)
1013  {
1014  if(r(i,i) == NumericTraits<T>::zero())
1015  return false; // r doesn't have full rank
1016  T sum = b(i, k);
1017  for(MultiArrayIndex j=i+1; j<m; ++j)
1018  sum -= r(i, j) * x(j, k);
1019  x(i, k) = sum / r(i, i);
1020  }
1021  }
1022  return true;
1023 }
1024 
1025  /** Solve a linear system with lower-triangular coefficient matrix.
1026 
1027  The square matrix \a l must be a lower-triangular coefficient matrix. If \a l
1028  doesn't have full rank the function fails and returns <tt>false</tt>,
1029  otherwise it returns <tt>true</tt>. The upper triangular part of matrix \a l will not be touched,
1030  so it doesn't need to contain zeros.
1031 
1032  The column vectors of matrix \a b are the right-hand sides of the equation (several equations
1033  with the same coefficients can thus be solved in one go). The result is returned
1034  in \a x, whose columns contain the solutions for the correspoinding
1035  columns of \a b. This implementation can be applied in-place, i.e. <tt>&b == &x</tt> is allowed.
1036  The following size requirements apply:
1037 
1038  \code
1039  rowCount(l) == columnCount(l);
1040  rowCount(l) == rowCount(b);
1041  columnCount(l) == rowCount(x);
1042  columnCount(b) == columnCount(x);
1043  \endcode
1044 
1045  <b>\#include</b> <<a href="linear__solve_8hxx-source.html">vigra/linear_solve.hxx</a>> or<br>
1046  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
1047  Namespaces: vigra and vigra::linalg
1048  */
1049 template <class T, class C1, class C2, class C3>
1052 {
1055  vigra_precondition(m == rowCount(l),
1056  "linearSolveLowerTriangular(): square coefficient matrix required.");
1057  vigra_precondition(m == rowCount(b) && m == rowCount(x) && n == columnCount(x),
1058  "linearSolveLowerTriangular(): matrix shape mismatch.");
1059 
1060  for(MultiArrayIndex k = 0; k < n; ++k)
1061  {
1062  for(MultiArrayIndex i=0; i<m; ++i)
1063  {
1064  if(l(i,i) == NumericTraits<T>::zero())
1065  return false; // l doesn't have full rank
1066  T sum = b(i, k);
1067  for(MultiArrayIndex j=0; j<i; ++j)
1068  sum -= l(i, j) * x(j, k);
1069  x(i, k) = sum / l(i, i);
1070  }
1071  }
1072  return true;
1073 }
1074 
1075 
1076  /** Solve a linear system when the Cholesky decomposition of the left hand side is given.
1077 
1078  The square matrix \a L must be a lower-triangular matrix resulting from Cholesky
1079  decomposition of some positive definite coefficient matrix.
1080 
1081  The column vectors of matrix \a b are the right-hand sides of the equation (several equations
1082  with the same matrix \a L can thus be solved in one go). The result is returned
1083  in \a x, whose columns contain the solutions for the correspoinding
1084  columns of \a b. This implementation can be applied in-place, i.e. <tt>&b == &x</tt> is allowed.
1085  The following size requirements apply:
1086 
1087  \code
1088  rowCount(L) == columnCount(L);
1089  rowCount(L) == rowCount(b);
1090  columnCount(L) == rowCount(x);
1091  columnCount(b) == columnCount(x);
1092  \endcode
1093 
1094  <b>\#include</b> <<a href="linear__solve_8hxx-source.html">vigra/linear_solve.hxx</a>> or<br>
1095  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
1096  Namespaces: vigra and vigra::linalg
1097  */
1098 template <class T, class C1, class C2, class C3>
1099 inline
1101 {
1102  /* Solve L * y = b */
1103  linearSolveLowerTriangular(L, b, x);
1104  /* Solve L^T * x = y */
1106 }
1107 
1108  /** Solve a linear system.
1109 
1110  \a A is the coefficient matrix, and the column vectors
1111  in \a b are the right-hand sides of the equation (so, several equations
1112  with the same coefficients can be solved in one go). The result is returned
1113  in \a res, whose columns contain the solutions for the corresponding
1114  columns of \a b. The number of columns of \a A must equal the number of rows of
1115  both \a b and \a res, and the number of columns of \a b and \a res must match.
1116 
1117  \a method must be one of the following:
1118  <DL>
1119  <DT>"Cholesky"<DD> Compute the solution by means of Cholesky decomposition. The
1120  coefficient matrix \a A must by symmetric positive definite. If
1121  this is not the case, the function returns <tt>false</tt>.
1122 
1123  <DT>"QR"<DD> (default) Compute the solution by means of QR decomposition. The
1124  coefficient matrix \a A can be square or rectangular. In the latter case,
1125  it must have more rows than columns, and the solution will be computed in the
1126  least squares sense. If \a A doesn't have full rank, the function
1127  returns <tt>false</tt>.
1128 
1129  <DT>"SVD"<DD> Compute the solution by means of singular value decomposition. The
1130  coefficient matrix \a A can be square or rectangular. In the latter case,
1131  it must have more rows than columns, and the solution will be computed in the
1132  least squares sense. If \a A doesn't have full rank, the function
1133  returns <tt>false</tt>.
1134 
1135  <DT>"NE"<DD> Compute the solution by means of the normal equations, i.e. by applying Cholesky
1136  decomposition to the equivalent problem <tt>A'*A*x = A'*b</tt>. This only makes sense
1137  when the equation is to be solved in the least squares sense, i.e. when \a A is a
1138  rectangular matrix with more rows than columns. If \a A doesn't have full column rank,
1139  the function returns <tt>false</tt>.
1140  </DL>
1141 
1142  This function can be applied in-place, i.e. <tt>&b == &res</tt> or <tt>&A == &res</tt> are allowed
1143  (provided they have the required shapes).
1144 
1145  The following size requirements apply:
1146 
1147  \code
1148  rowCount(r) == rowCount(b);
1149  columnCount(r) == rowCount(x);
1150  columnCount(b) == columnCount(x);
1151  \endcode
1152 
1153  <b>\#include</b> <<a href="linear__solve_8hxx-source.html">vigra/linear_solve.hxx</a>> or<br>
1154  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
1155  Namespaces: vigra and vigra::linalg
1156  */
1157 template <class T, class C1, class C2, class C3>
1159  MultiArrayView<2, T, C3> & res, std::string method = "QR")
1160 {
1161  typedef typename Matrix<T>::difference_type Shape;
1162  typedef typename Matrix<T>::view_type SubMatrix;
1163 
1164  const MultiArrayIndex n = columnCount(A);
1165  const MultiArrayIndex m = rowCount(A);
1166 
1167  vigra_precondition(n <= m,
1168  "linearSolve(): Coefficient matrix A must have at least as many rows as columns.");
1169  vigra_precondition(n == rowCount(res) &&
1170  m == rowCount(b) && columnCount(b) == columnCount(res),
1171  "linearSolve(): matrix shape mismatch.");
1172 
1173  for(unsigned int k=0; k<method.size(); ++k)
1174  method[k] = (std::string::value_type)tolower(method[k]);
1175 
1176  if(method == "cholesky")
1177  {
1178  vigra_precondition(columnCount(A) == rowCount(A),
1179  "linearSolve(): Cholesky method requires square coefficient matrix.");
1180  Matrix<T> L(A.shape());
1181  if(!choleskyDecomposition(A, L))
1182  return false; // false if A wasn't symmetric positive definite
1183  choleskySolve(L, b, res);
1184  }
1185  else if(method == "qr")
1186  {
1187  return (MultiArrayIndex)linearSolveQR(A, b, res) == n;
1188  }
1189  else if(method == "ne")
1190  {
1191  return linearSolve(transpose(A)*A, transpose(A)*b, res, "Cholesky");
1192  }
1193  else if(method == "svd")
1194  {
1195  MultiArrayIndex rhsCount = columnCount(b);
1196  Matrix<T> u(A.shape()), s(n, 1), v(n, n);
1197 
1199 
1200  Matrix<T> t = transpose(u)*b;
1201  for(MultiArrayIndex l=0; l<rhsCount; ++l)
1202  {
1203  for(MultiArrayIndex k=0; k<rank; ++k)
1204  t(k,l) /= s(k,0);
1205  for(MultiArrayIndex k=rank; k<n; ++k)
1206  t(k,l) = NumericTraits<T>::zero();
1207  }
1208  res = v*t;
1209 
1210  return (rank == n);
1211  }
1212  else
1213  {
1214  vigra_precondition(false, "linearSolve(): Unknown solution method.");
1215  }
1216  return true;
1217 }
1218 
1219 //@}
1220 
1221 } // namespace linalg
1222 
1223 using linalg::inverse;
1224 using linalg::determinant;
1226 using linalg::linearSolve;
1227 using linalg::choleskySolve;
1232 
1233 } // namespace vigra
1234 
1235 
1236 #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.7.1 (Thu Jun 14 2012)