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

linear_solve.hxx VIGRA

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>static_cast<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>static_cast<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), static_cast<int>(n)),
394  columnVector(z, Shape(0,0), static_cast<int>(n)));
395  // use atan2 as it is robust against overflow/underflow
396  T t = 0.5*std::atan2(T(-2.0*yv), T(squaredNorm(gamma / v) + squaredNorm(yv) - 1.0)),
397  s = std::sin(t),
398  c = std::cos(t);
399  columnVector(z, Shape(0,0), static_cast<int>(n)) *= c;
400  z(n,0) = (s - c*yv) / gamma;
401  v *= norm(gamma) / hypot(c*gamma, v*(s - c*yv));
402 }
403 
404 // QR algorithm with optional column pivoting
405 template <class T, class C1, class C2, class C3>
406 unsigned int
407 qrTransformToTriangularImpl(MultiArrayView<2, T, C1> & r, MultiArrayView<2, T, C2> & rhs, MultiArrayView<2, T, C3> & householder,
408  ArrayVector<MultiArrayIndex> & permutation, double epsilon)
409 {
410  typedef typename Matrix<T>::difference_type Shape;
411  typedef typename NormTraits<MultiArrayView<2, T, C1> >::NormType NormType;
412  typedef typename NormTraits<MultiArrayView<2, T, C1> >::SquaredNormType SNType;
413 
414  const MultiArrayIndex m = rowCount(r);
415  const MultiArrayIndex n = columnCount(r);
416  const MultiArrayIndex maxRank = std::min(m, n);
417 
418  vigra_precondition(m >= n,
419  "qrTransformToTriangularImpl(): Coefficient matrix with at least as many rows as columns required.");
420 
421  const MultiArrayIndex rhsCount = columnCount(rhs);
422  bool transformRHS = rhsCount > 0;
423  vigra_precondition(!transformRHS || m == rowCount(rhs),
424  "qrTransformToTriangularImpl(): RHS matrix shape mismatch.");
425 
426  bool storeHouseholderSteps = columnCount(householder) > 0;
427  vigra_precondition(!storeHouseholderSteps || r.shape() == householder.shape(),
428  "qrTransformToTriangularImpl(): Householder matrix shape mismatch.");
429 
430  bool pivoting = permutation.size() > 0;
431  vigra_precondition(!pivoting || n == static_cast<MultiArrayIndex>(permutation.size()),
432  "qrTransformToTriangularImpl(): Permutation array size mismatch.");
433 
434  if(n == 0)
435  return 0; // trivial solution
436 
437  Matrix<SNType> columnSquaredNorms;
438  if(pivoting)
439  {
440  columnSquaredNorms.reshape(Shape(1,n));
441  for(MultiArrayIndex k=0; k<n; ++k)
442  columnSquaredNorms[k] = squaredNorm(columnVector(r, k));
443 
444  int pivot = argMax(columnSquaredNorms);
445  if(pivot != 0)
446  {
447  columnVector(r, 0).swapData(columnVector(r, pivot));
448  std::swap(columnSquaredNorms[0], columnSquaredNorms[pivot]);
449  std::swap(permutation[0], permutation[pivot]);
450  }
451  }
452 
453  qrHouseholderStepImpl(0, r, rhs, householder);
454 
455  MultiArrayIndex rank = 1;
456  NormType maxApproxSingularValue = norm(r(0,0)),
457  minApproxSingularValue = maxApproxSingularValue;
458 
459  double tolerance = (epsilon == 0.0)
460  ? m*maxApproxSingularValue*NumericTraits<T>::epsilon()
461  : epsilon;
462 
463  bool simpleSingularValueApproximation = (n < 4);
464  Matrix<T> zmax, zmin;
465  if(minApproxSingularValue <= tolerance)
466  {
467  rank = 0;
468  pivoting = false;
469  simpleSingularValueApproximation = true;
470  }
471  if(!simpleSingularValueApproximation)
472  {
473  zmax.reshape(Shape(m,1));
474  zmin.reshape(Shape(m,1));
475  zmax(0,0) = r(0,0);
476  zmin(0,0) = 1.0 / r(0,0);
477  }
478 
479  for(MultiArrayIndex k=1; k<maxRank; ++k)
480  {
481  if(pivoting)
482  {
483  for(MultiArrayIndex l=k; l<n; ++l)
484  columnSquaredNorms[l] -= squaredNorm(r(k, l));
485  int pivot = k + argMax(rowVector(columnSquaredNorms, Shape(0,k), n));
486  if(pivot != static_cast<int>(k))
487  {
488  columnVector(r, k).swapData(columnVector(r, pivot));
489  std::swap(columnSquaredNorms[k], columnSquaredNorms[pivot]);
490  std::swap(permutation[k], permutation[pivot]);
491  }
492  }
493 
494  qrHouseholderStepImpl(k, r, rhs, householder);
495 
496  if(simpleSingularValueApproximation)
497  {
498  NormType nv = norm(r(k,k));
499  maxApproxSingularValue = std::max(nv, maxApproxSingularValue);
500  minApproxSingularValue = std::min(nv, minApproxSingularValue);
501  }
502  else
503  {
504  incrementalMaxSingularValueApproximation(columnVector(r, Shape(0,k),k+1), zmax, maxApproxSingularValue);
505  incrementalMinSingularValueApproximation(columnVector(r, Shape(0,k),k+1), zmin, minApproxSingularValue, tolerance);
506  }
507 
508 #if 0
509  Matrix<T> u(k+1,k+1), s(k+1, 1), v(k+1,k+1);
510  singularValueDecomposition(r.subarray(Shape(0,0), Shape(k+1,k+1)), u, s, v);
511  std::cerr << "estimate, svd " << k << ": " << minApproxSingularValue << " " << s(k,0) << "\n";
512 #endif
513 
514  if(epsilon == 0.0)
515  tolerance = m*maxApproxSingularValue*NumericTraits<T>::epsilon();
516 
517  if(minApproxSingularValue > tolerance)
518  ++rank;
519  else
520  pivoting = false; // matrix doesn't have full rank, triangulize the rest without pivoting
521  }
522  return static_cast<unsigned int>(rank);
523 }
524 
525 template <class T, class C1, class C2>
526 unsigned int
527 qrTransformToUpperTriangular(MultiArrayView<2, T, C1> & r, MultiArrayView<2, T, C2> & rhs,
528  ArrayVector<MultiArrayIndex> & permutation, double epsilon = 0.0)
529 {
530  Matrix<T> dontStoreHouseholderVectors; // intentionally empty
531  return qrTransformToTriangularImpl(r, rhs, dontStoreHouseholderVectors, permutation, epsilon);
532 }
533 
534 // QR algorithm with optional row pivoting
535 template <class T, class C1, class C2, class C3>
536 unsigned int
537 qrTransformToLowerTriangular(MultiArrayView<2, T, C1> & r, MultiArrayView<2, T, C2> & rhs, MultiArrayView<2, T, C3> & householderMatrix,
538  double epsilon = 0.0)
539 {
540  ArrayVector<MultiArrayIndex> permutation(static_cast<unsigned int>(rowCount(rhs)));
541  for(MultiArrayIndex k=0; k<static_cast<MultiArrayIndex>(permutation.size()); ++k)
542  permutation[k] = k;
543  Matrix<T> dontTransformRHS; // intentionally empty
544  MultiArrayView<2, T, StridedArrayTag> rt = transpose(r),
545  ht = transpose(householderMatrix);
546  unsigned int rank = qrTransformToTriangularImpl(rt, dontTransformRHS, ht, permutation, epsilon);
547 
548  // apply row permutation to RHS
549  Matrix<T> tempRHS(rhs);
550  for(MultiArrayIndex k=0; k<static_cast<MultiArrayIndex>(permutation.size()); ++k)
551  rowVector(rhs, k) = rowVector(tempRHS, permutation[k]);
552  return rank;
553 }
554 
555 // QR algorithm without column pivoting
556 template <class T, class C1, class C2>
557 inline bool
558 qrTransformToUpperTriangular(MultiArrayView<2, T, C1> & r, MultiArrayView<2, T, C2> & rhs,
559  double epsilon = 0.0)
560 {
561  ArrayVector<MultiArrayIndex> noPivoting; // intentionally empty
562 
563  return (qrTransformToUpperTriangular(r, rhs, noPivoting, epsilon) ==
564  static_cast<unsigned int>(columnCount(r)));
565 }
566 
567 // QR algorithm without row pivoting
568 template <class T, class C1, class C2>
569 inline bool
570 qrTransformToLowerTriangular(MultiArrayView<2, T, C1> & r, MultiArrayView<2, T, C2> & householder,
571  double epsilon = 0.0)
572 {
573  Matrix<T> noPivoting; // intentionally empty
574 
575  return (qrTransformToLowerTriangular(r, noPivoting, householder, epsilon) ==
576  static_cast<unsigned int>(rowCount(r)));
577 }
578 
579 // restore ordering of result vector elements after QR solution with column pivoting
580 template <class T, class C1, class C2, class Permutation>
581 void inverseRowPermutation(MultiArrayView<2, T, C1> &permuted, MultiArrayView<2, T, C2> &res,
582  Permutation const & permutation)
583 {
584  for(MultiArrayIndex k=0; k<columnCount(permuted); ++k)
585  for(MultiArrayIndex l=0; l<rowCount(permuted); ++l)
586  res(permutation[l], k) = permuted(l,k);
587 }
588 
589 template <class T, class C1, class C2>
590 void applyHouseholderColumnReflections(MultiArrayView<2, T, C1> const &householder, MultiArrayView<2, T, C2> &res)
591 {
592  typedef typename Matrix<T>::difference_type Shape;
593  MultiArrayIndex n = rowCount(householder);
594  MultiArrayIndex m = columnCount(householder);
595  MultiArrayIndex rhsCount = columnCount(res);
596 
597  for(int k = m-1; k >= 0; --k)
598  {
599  MultiArrayView<2, T, C1> u = columnVector(householder, Shape(k,k), n);
600  for(MultiArrayIndex l=0; l<rhsCount; ++l)
601  columnVector(res, Shape(k,l), n) -= dot(columnVector(res, Shape(k,l), n), u) * u;
602  }
603 }
604 
605 } // namespace detail
606 
607 template <class T, class C1, class C2, class C3>
608 unsigned int
609 linearSolveQRReplace(MultiArrayView<2, T, C1> &A, MultiArrayView<2, T, C2> &b,
610  MultiArrayView<2, T, C3> & res,
611  double epsilon = 0.0)
612 {
613  typedef typename Matrix<T>::difference_type Shape;
614 
616  MultiArrayIndex m = rowCount(A);
617  MultiArrayIndex rhsCount = columnCount(res);
618  MultiArrayIndex rank = std::min(m,n);
619  Shape ul(MultiArrayIndex(0), MultiArrayIndex(0));
620 
621 
622  vigra_precondition(rhsCount == columnCount(b),
623  "linearSolveQR(): RHS and solution must have the same number of columns.");
624  vigra_precondition(m == rowCount(b),
625  "linearSolveQR(): Coefficient matrix and RHS must have the same number of rows.");
626  vigra_precondition(n == rowCount(res),
627  "linearSolveQR(): Mismatch between column count of coefficient matrix and row count of solution.");
628  vigra_precondition(epsilon >= 0.0,
629  "linearSolveQR(): 'epsilon' must be non-negative.");
630 
631  if(m < n)
632  {
633  // minimum norm solution of underdetermined system
634  Matrix<T> householderMatrix(n, m);
635  MultiArrayView<2, T, StridedArrayTag> ht = transpose(householderMatrix);
636  rank = static_cast<MultiArrayIndex>(detail::qrTransformToLowerTriangular(A, b, ht, epsilon));
637  res.subarray(Shape(rank,0), Shape(n, rhsCount)).init(NumericTraits<T>::zero());
638  if(rank < m)
639  {
640  // system is also rank-deficient => compute minimum norm least squares solution
641  MultiArrayView<2, T, C1> Asub = A.subarray(ul, Shape(m,rank));
642  detail::qrTransformToUpperTriangular(Asub, b, epsilon);
643  linearSolveUpperTriangular(A.subarray(ul, Shape(rank,rank)),
644  b.subarray(ul, Shape(rank,rhsCount)),
645  res.subarray(ul, Shape(rank, rhsCount)));
646  }
647  else
648  {
649  // system has full rank => compute minimum norm solution
650  linearSolveLowerTriangular(A.subarray(ul, Shape(rank,rank)),
651  b.subarray(ul, Shape(rank, rhsCount)),
652  res.subarray(ul, Shape(rank, rhsCount)));
653  }
654  detail::applyHouseholderColumnReflections(householderMatrix.subarray(ul, Shape(n, rank)), res);
655  }
656  else
657  {
658  // solution of well-determined or overdetermined system
659  ArrayVector<MultiArrayIndex> permutation(static_cast<unsigned int>(n));
660  for(MultiArrayIndex k=0; k<n; ++k)
661  permutation[k] = k;
662 
663  rank = static_cast<MultiArrayIndex>(detail::qrTransformToUpperTriangular(A, b, permutation, epsilon));
664 
665  Matrix<T> permutedSolution(n, rhsCount);
666  if(rank < n)
667  {
668  // system is rank-deficient => compute minimum norm solution
669  Matrix<T> householderMatrix(n, rank);
670  MultiArrayView<2, T, StridedArrayTag> ht = transpose(householderMatrix);
671  MultiArrayView<2, T, C1> Asub = A.subarray(ul, Shape(rank,n));
672  detail::qrTransformToLowerTriangular(Asub, ht, epsilon);
673  linearSolveLowerTriangular(A.subarray(ul, Shape(rank,rank)),
674  b.subarray(ul, Shape(rank, rhsCount)),
675  permutedSolution.subarray(ul, Shape(rank, rhsCount)));
676  detail::applyHouseholderColumnReflections(householderMatrix, permutedSolution);
677  }
678  else
679  {
680  // system has full rank => compute exact or least squares solution
681  linearSolveUpperTriangular(A.subarray(ul, Shape(rank,rank)),
682  b.subarray(ul, Shape(rank,rhsCount)),
683  permutedSolution);
684  }
685  detail::inverseRowPermutation(permutedSolution, res, permutation);
686  }
687  return static_cast<unsigned int>(rank);
688 }
689 
690 template <class T, class C1, class C2, class C3>
691 unsigned int linearSolveQR(MultiArrayView<2, T, C1> const & A, MultiArrayView<2, T, C2> const & b,
692  MultiArrayView<2, T, C3> & res)
693 {
694  Matrix<T> r(A), rhs(b);
695  return linearSolveQRReplace(r, rhs, res);
696 }
697 
698 /** \defgroup MatrixAlgebra Advanced Matrix Algebra
699 
700  \brief Solution of linear systems, eigen systems, linear least squares etc.
701 
702  \ingroup LinearAlgebraModule
703  */
704 //@{
705  /** Create the inverse or pseudo-inverse of matrix \a v.
706 
707  If the matrix \a v is square, \a res must have the same shape and will contain the
708  inverse of \a v. If \a v is rectangular, \a res must have the transposed shape
709  of \a v. The inverse is then computed in the least-squares
710  sense, i.e. \a res will be the pseudo-inverse (Moore-Penrose inverse).
711  The function returns <tt>true</tt> upon success, and <tt>false</tt> if \a v
712  is not invertible (has not full rank). The inverse is computed by means of QR
713  decomposition. This function can be applied in-place.
714 
715  <b>\#include</b> <vigra/linear_solve.hxx> or<br>
716  <b>\#include</b> <vigra/linear_algebra.hxx><br>
717  Namespaces: vigra and vigra::linalg
718  */
719 template <class T, class C1, class C2>
721 {
722  typedef typename MultiArrayShape<2>::type Shape;
723 
724  const MultiArrayIndex n = columnCount(v);
725  const MultiArrayIndex m = rowCount(v);
726  vigra_precondition(n == rowCount(res) && m == columnCount(res),
727  "inverse(): shape of output matrix must be the transpose of the input matrix' shape.");
728 
729  if(m < n)
730  {
732  Matrix<T> r(vt.shape()), q(n, n);
733  if(!qrDecomposition(vt, q, r))
734  return false; // a didn't have full rank
735  linearSolveUpperTriangular(r.subarray(Shape(0,0), Shape(m,m)),
736  transpose(q).subarray(Shape(0,0), Shape(m,n)),
737  transpose(res));
738  }
739  else
740  {
741  Matrix<T> r(v.shape()), q(m, m);
742  if(!qrDecomposition(v, q, r))
743  return false; // a didn't have full rank
744  linearSolveUpperTriangular(r.subarray(Shape(0,0), Shape(n,n)),
745  transpose(q).subarray(Shape(0,0), Shape(n,m)),
746  res);
747  }
748  return true;
749 }
750 
751  /** Create the inverse or pseudo-inverse of matrix \a v.
752 
753  The result is returned as a temporary matrix. If the matrix \a v is square,
754  the result will have the same shape and contains the inverse of \a v.
755  If \a v is rectangular, the result will have the transposed shape of \a v.
756  The inverse is then computed in the least-squares
757  sense, i.e. \a res will be the pseudo-inverse (Moore-Penrose inverse).
758  The inverse is computed by means of QR decomposition. If \a v
759  is not invertible, <tt>vigra::PreconditionViolation</tt> exception is thrown.
760  Usage:
761 
762  \code
763  vigra::Matrix<double> v(n, n);
764  v = ...;
765 
766  vigra::Matrix<double> m = inverse(v);
767  \endcode
768 
769  <b>\#include</b> <vigra/linear_solve.hxx> or<br>
770  <b>\#include</b> <vigra/linear_algebra.hxx><br>
771  Namespaces: vigra and vigra::linalg
772  */
773 template <class T, class C>
774 TemporaryMatrix<T> inverse(const MultiArrayView<2, T, C> &v)
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 template <class T>
783 TemporaryMatrix<T> inverse(const TemporaryMatrix<T> &v)
784 {
785  if(columnCount(v) == rowCount(v))
786  {
787  vigra_precondition(inverse(v, const_cast<TemporaryMatrix<T> &>(v)),
788  "inverse(): matrix is not invertible.");
789  return v;
790  }
791  else
792  {
793  TemporaryMatrix<T> ret(columnCount(v), rowCount(v)); // transpose shape
794  vigra_precondition(inverse(v, ret),
795  "inverse(): matrix is not invertible.");
796  return ret;
797  }
798 }
799 
800  /** Compute the determinant of a square matrix.
801 
802  \a method must be one of the following:
803  <DL>
804  <DT>"Cholesky"<DD> Compute the solution by means of Cholesky decomposition. This
805  method is faster than "LU", but requires the matrix \a a
806  to be symmetric positive definite. If this is
807  not the case, a <tt>ContractViolation</tt> exception is thrown.
808 
809  <DT>"LU"<DD> (default) Compute the solution by means of LU decomposition.
810  </DL>
811 
812  <b>\#include</b> <vigra/linear_solve.hxx> or<br>
813  <b>\#include</b> <vigra/linear_algebra.hxx><br>
814  Namespaces: vigra and vigra::linalg
815  */
816 template <class T, class C1>
817 T determinant(MultiArrayView<2, T, C1> const & a, std::string method = "LU")
818 {
820  vigra_precondition(rowCount(a) == n,
821  "determinant(): Square matrix required.");
822 
823  method = tolower(method);
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 (static_cast<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  MultiArrayIndex m = rowCount(r);
1019  MultiArrayIndex rhsCount = columnCount(b);
1020  vigra_precondition(m == columnCount(r),
1021  "linearSolveUpperTriangular(): square coefficient matrix required.");
1022  vigra_precondition(m == rowCount(b) && m == rowCount(x) && rhsCount == columnCount(x),
1023  "linearSolveUpperTriangular(): matrix shape mismatch.");
1024 
1025  for(MultiArrayIndex k = 0; k < rhsCount; ++k)
1026  {
1027  for(int i=m-1; i>=0; --i)
1028  {
1029  if(r(i,i) == NumericTraits<T>::zero())
1030  return false; // r doesn't have full rank
1031  T sum = b(i, k);
1032  for(MultiArrayIndex j=i+1; j<m; ++j)
1033  sum -= r(i, j) * x(j, k);
1034  x(i, k) = sum / r(i, i);
1035  }
1036  }
1037  return true;
1038 }
1039 
1040  /** Solve a linear system with lower-triangular coefficient matrix.
1041 
1042  The square matrix \a l must be a lower-triangular coefficient matrix. If \a l
1043  doesn't have full rank the function fails and returns <tt>false</tt>,
1044  otherwise it returns <tt>true</tt>. The upper triangular part of matrix \a l will not be touched,
1045  so it doesn't need to contain zeros.
1046 
1047  The column vectors of matrix \a b are the right-hand sides of the equation (several equations
1048  with the same coefficients can thus be solved in one go). The result is returned
1049  in \a x, whose columns contain the solutions for the corresponding
1050  columns of \a b. This implementation can be applied in-place, i.e. <tt>&b == &x</tt> is allowed.
1051  The following size requirements apply:
1052 
1053  \code
1054  rowCount(l) == columnCount(l);
1055  rowCount(l) == rowCount(b);
1056  columnCount(l) == rowCount(x);
1057  columnCount(b) == columnCount(x);
1058  \endcode
1059 
1060  <b>\#include</b> <vigra/linear_solve.hxx> or<br>
1061  <b>\#include</b> <vigra/linear_algebra.hxx><br>
1062  Namespaces: vigra and vigra::linalg
1063  */
1064 template <class T, class C1, class C2, class C3>
1067 {
1070  vigra_precondition(m == rowCount(l),
1071  "linearSolveLowerTriangular(): square coefficient matrix required.");
1072  vigra_precondition(m == rowCount(b) && m == rowCount(x) && n == columnCount(x),
1073  "linearSolveLowerTriangular(): matrix shape mismatch.");
1074 
1075  for(MultiArrayIndex k = 0; k < n; ++k)
1076  {
1077  for(MultiArrayIndex i=0; i<m; ++i)
1078  {
1079  if(l(i,i) == NumericTraits<T>::zero())
1080  return false; // l doesn't have full rank
1081  T sum = b(i, k);
1082  for(MultiArrayIndex j=0; j<i; ++j)
1083  sum -= l(i, j) * x(j, k);
1084  x(i, k) = sum / l(i, i);
1085  }
1086  }
1087  return true;
1088 }
1089 
1090 
1091  /** Solve a linear system when the Cholesky decomposition of the left hand side is given.
1092 
1093  The square matrix \a L must be a lower-triangular matrix resulting from Cholesky
1094  decomposition of some positive definite coefficient matrix.
1095 
1096  The column vectors of matrix \a b are the right-hand sides of the equation (several equations
1097  with the same matrix \a L can thus be solved in one go). The result is returned
1098  in \a x, whose columns contain the solutions for the corresponding
1099  columns of \a b. This implementation can be applied in-place, i.e. <tt>&b == &x</tt> is allowed.
1100  The following size requirements apply:
1101 
1102  \code
1103  rowCount(L) == columnCount(L);
1104  rowCount(L) == rowCount(b);
1105  columnCount(L) == rowCount(x);
1106  columnCount(b) == columnCount(x);
1107  \endcode
1108 
1109  <b>\#include</b> <vigra/linear_solve.hxx> or<br>
1110  <b>\#include</b> <vigra/linear_algebra.hxx><br>
1111  Namespaces: vigra and vigra::linalg
1112  */
1113 template <class T, class C1, class C2, class C3>
1114 inline
1116 {
1117  /* Solve L * y = b */
1118  linearSolveLowerTriangular(L, b, x);
1119  /* Solve L^T * x = y */
1121 }
1122 
1123  /** Solve a linear system.
1124 
1125  <b> Declarations:</b>
1126 
1127  \code
1128  // use MultiArrayViews for input and output
1129  template <class T, class C1, class C2, class C3>
1130  bool linearSolve(MultiArrayView<2, T, C1> const & A,
1131  MultiArrayView<2, T, C2> const & b,
1132  MultiArrayView<2, T, C3> res,
1133  std::string method = "QR");
1134 
1135  // use TinyVector for RHS and result
1136  template <class T, class C1, int N>
1137  bool linearSolve(MultiArrayView<2, T, C1> const & A,
1138  TinyVector<T, N> const & b,
1139  TinyVector<T, N> & res,
1140  std::string method = "QR");
1141  \endcode
1142 
1143  \a A is the coefficient matrix, and the column vectors
1144  in \a b are the right-hand sides of the equation (so, several equations
1145  with the same coefficients can be solved in one go). The result is returned
1146  in \a res, whose columns contain the solutions for the corresponding
1147  columns of \a b. The number of columns of \a A must equal the number of rows of
1148  both \a b and \a res, and the number of columns of \a b and \a res must match.
1149  If right-hand-side and result are specified as TinyVector, the number of columns
1150  of \a A must equal N.
1151 
1152  \a method must be one of the following:
1153  <DL>
1154  <DT>"Cholesky"<DD> Compute the solution by means of Cholesky decomposition. The
1155  coefficient matrix \a A must by symmetric positive definite. If
1156  this is not the case, the function returns <tt>false</tt>.
1157 
1158  <DT>"QR"<DD> (default) Compute the solution by means of QR decomposition. The
1159  coefficient matrix \a A can be square or rectangular. In the latter case,
1160  it must have more rows than columns, and the solution will be computed in the
1161  least squares sense. If \a A doesn't have full rank, the function
1162  returns <tt>false</tt>.
1163 
1164  <DT>"SVD"<DD> Compute the solution by means of singular value decomposition. The
1165  coefficient matrix \a A can be square or rectangular. In the latter case,
1166  it must have more rows than columns, and the solution will be computed in the
1167  least squares sense. If \a A doesn't have full rank, the function
1168  returns <tt>false</tt>.
1169 
1170  <DT>"NE"<DD> Compute the solution by means of the normal equations, i.e. by applying Cholesky
1171  decomposition to the equivalent problem <tt>A'*A*x = A'*b</tt>. This only makes sense
1172  when the equation is to be solved in the least squares sense, i.e. when \a A is a
1173  rectangular matrix with more rows than columns. If \a A doesn't have full column rank,
1174  the function returns <tt>false</tt>.
1175  </DL>
1176 
1177  This function can be applied in-place, i.e. <tt>&b == &res</tt> or <tt>&A == &res</tt> are allowed
1178  (provided they have the required shapes).
1179 
1180  The following size requirements apply:
1181 
1182  \code
1183  rowCount(r) == rowCount(b);
1184  columnCount(r) == rowCount(x);
1185  columnCount(b) == columnCount(x);
1186  \endcode
1187 
1188  <b>\#include</b> <vigra/linear_solve.hxx> or<br>
1189  <b>\#include</b> <vigra/linear_algebra.hxx><br>
1190  Namespaces: vigra and vigra::linalg
1191  */
1192 doxygen_overloaded_function(template <...> bool linearSolve)
1193 
1194 template <class T, class C1, class C2, class C3>
1195 bool linearSolve(MultiArrayView<2, T, C1> const & A,
1196  MultiArrayView<2, T, C2> const & b,
1198  std::string method = "QR")
1199 {
1200  const MultiArrayIndex n = columnCount(A);
1201  const MultiArrayIndex m = rowCount(A);
1202 
1203  vigra_precondition(n <= m,
1204  "linearSolve(): Coefficient matrix A must have at least as many rows as columns.");
1205  vigra_precondition(n == rowCount(res) &&
1206  m == rowCount(b) && columnCount(b) == columnCount(res),
1207  "linearSolve(): matrix shape mismatch.");
1208 
1209  method = tolower(method);
1210  if(method == "cholesky")
1211  {
1212  vigra_precondition(columnCount(A) == rowCount(A),
1213  "linearSolve(): Cholesky method requires square coefficient matrix.");
1214  Matrix<T> L(A.shape());
1215  if(!choleskyDecomposition(A, L))
1216  return false; // false if A wasn't symmetric positive definite
1217  choleskySolve(L, b, res);
1218  }
1219  else if(method == "qr")
1220  {
1221  return static_cast<MultiArrayIndex>(linearSolveQR(A, b, res)) == n;
1222  }
1223  else if(method == "ne")
1224  {
1225  return linearSolve(transpose(A)*A, transpose(A)*b, res, "Cholesky");
1226  }
1227  else if(method == "svd")
1228  {
1229  MultiArrayIndex rhsCount = columnCount(b);
1230  Matrix<T> u(A.shape()), s(n, 1), v(n, n);
1231 
1232  MultiArrayIndex rank = static_cast<MultiArrayIndex>(singularValueDecomposition(A, u, s, v));
1233 
1234  Matrix<T> t = transpose(u)*b;
1235  for(MultiArrayIndex l=0; l<rhsCount; ++l)
1236  {
1237  for(MultiArrayIndex k=0; k<rank; ++k)
1238  t(k,l) /= s(k,0);
1239  for(MultiArrayIndex k=rank; k<n; ++k)
1240  t(k,l) = NumericTraits<T>::zero();
1241  }
1242  res = v*t;
1243 
1244  return (rank == n);
1245  }
1246  else
1247  {
1248  vigra_precondition(false, "linearSolve(): Unknown solution method.");
1249  }
1250  return true;
1251 }
1252 
1253 template <class T, class C1, int N>
1254 bool linearSolve(MultiArrayView<2, T, C1> const & A,
1255  TinyVector<T, N> const & b,
1256  TinyVector<T, N> & res,
1257  std::string method = "QR")
1258 {
1259  Shape2 shape(N, 1);
1260  return linearSolve(A, MultiArrayView<2, T>(shape, b.data()), MultiArrayView<2, T>(shape, res.data()), method);
1261 }
1262 
1263 //@}
1264 
1265 } // namespace linalg
1266 
1267 using linalg::inverse;
1268 using linalg::determinant;
1270 using linalg::linearSolve;
1271 using linalg::choleskySolve;
1276 
1277 } // namespace vigra
1278 
1279 
1280 #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.10.0 (Thu Jan 8 2015)