[ 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>(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  method = tolower(method);
823 
824  if(n == 1)
825  return a(0,0);
826  if(n == 2)
827  return a(0,0)*a(1,1) - a(0,1)*a(1,0);
828  if(method == "lu")
829  {
830  return detail::determinantByLUDecomposition(a);
831  }
832  else if(method == "cholesky")
833  {
834  Matrix<T> L(a.shape());
835  vigra_precondition(choleskyDecomposition(a, L),
836  "determinant(): Cholesky method requires symmetric positive definite matrix.");
837  T det = L(0,0);
838  for(MultiArrayIndex k=1; k<n; ++k)
839  det *= L(k,k);
840  return sq(det);
841  }
842  else
843  {
844  vigra_precondition(false, "determinant(): Unknown solution method.");
845  }
846  return T();
847 }
848 
849  /** Compute the logarithm of the determinant of a symmetric positive definite matrix.
850 
851  This is useful to avoid multiplication of very large numbers in big matrices.
852  It is implemented by means of Cholesky decomposition.
853 
854  <b>\#include</b> <vigra/linear_solve.hxx> or<br>
855  <b>\#include</b> <vigra/linear_algebra.hxx><br>
856  Namespaces: vigra and vigra::linalg
857  */
858 template <class T, class C1>
860 {
862  vigra_precondition(rowCount(a) == n,
863  "logDeterminant(): Square matrix required.");
864  if(n == 1)
865  {
866  vigra_precondition(a(0,0) > 0.0,
867  "logDeterminant(): Matrix not positive definite.");
868  return std::log(a(0,0));
869  }
870  if(n == 2)
871  {
872  T det = a(0,0)*a(1,1) - a(0,1)*a(1,0);
873  vigra_precondition(det > 0.0,
874  "logDeterminant(): Matrix not positive definite.");
875  return std::log(det);
876  }
877  else
878  {
879  Matrix<T> L(a.shape());
880  vigra_precondition(choleskyDecomposition(a, L),
881  "logDeterminant(): Matrix not positive definite.");
882  T logdet = std::log(L(0,0));
883  for(MultiArrayIndex k=1; k<n; ++k)
884  logdet += std::log(L(k,k)); // L(k,k) is guaranteed to be positive
885  return 2.0*logdet;
886  }
887 }
888 
889  /** Cholesky decomposition.
890 
891  \a A must be a symmetric positive definite matrix, and \a L will be a lower
892  triangular matrix, such that (up to round-off errors):
893 
894  \code
895  A == L * transpose(L);
896  \endcode
897 
898  This implementation cannot be applied in-place, i.e. <tt>&L == &A</tt> is an error.
899  If \a A is not symmetric, a <tt>ContractViolation</tt> exception is thrown. If it
900  is not positive definite, the function returns <tt>false</tt>.
901 
902  <b>\#include</b> <vigra/linear_solve.hxx> or<br>
903  <b>\#include</b> <vigra/linear_algebra.hxx><br>
904  Namespaces: vigra and vigra::linalg
905  */
906 template <class T, class C1, class C2>
909 {
911  vigra_precondition(rowCount(A) == n,
912  "choleskyDecomposition(): Input matrix must be square.");
913  vigra_precondition(n == columnCount(L) && n == rowCount(L),
914  "choleskyDecomposition(): Output matrix must have same shape as input matrix.");
915  vigra_precondition(isSymmetric(A),
916  "choleskyDecomposition(): Input matrix must be symmetric.");
917 
918  for (MultiArrayIndex j = 0; j < n; ++j)
919  {
920  T d(0.0);
921  for (MultiArrayIndex k = 0; k < j; ++k)
922  {
923  T s(0.0);
924  for (MultiArrayIndex i = 0; i < k; ++i)
925  {
926  s += L(k, i)*L(j, i);
927  }
928  L(j, k) = s = (A(j, k) - s)/L(k, k);
929  d = d + s*s;
930  }
931  d = A(j, j) - d;
932  if(d <= 0.0)
933  return false; // A is not positive definite
934  L(j, j) = std::sqrt(d);
935  for (MultiArrayIndex k = j+1; k < n; ++k)
936  {
937  L(j, k) = 0.0;
938  }
939  }
940  return true;
941 }
942 
943  /** QR decomposition.
944 
945  \a a contains the original matrix, results are returned in \a q and \a r, where
946  \a q is a orthogonal matrix, and \a r is an upper triangular matrix, such that
947  (up to round-off errors):
948 
949  \code
950  a == q * r;
951  \endcode
952 
953  If \a a doesn't have full rank, the function returns <tt>false</tt>.
954  The decomposition is computed by householder transformations. It can be applied in-place,
955  i.e. <tt>&a == &q</tt> or <tt>&a == &r</tt> are allowed.
956 
957  <b>\#include</b> <vigra/linear_solve.hxx> or<br>
958  <b>\#include</b> <vigra/linear_algebra.hxx><br>
959  Namespaces: vigra and vigra::linalg
960  */
961 template <class T, class C1, class C2, class C3>
964  double epsilon = 0.0)
965 {
966  const MultiArrayIndex m = rowCount(a);
967  const MultiArrayIndex n = columnCount(a);
968  vigra_precondition(n == columnCount(r) && m == rowCount(r) &&
969  m == columnCount(q) && m == rowCount(q),
970  "qrDecomposition(): Matrix shape mismatch.");
971 
972  q = identityMatrix<T>(m);
974  r = a;
975  ArrayVector<MultiArrayIndex> noPivoting; // intentionally empty
976  return ((MultiArrayIndex)detail::qrTransformToUpperTriangular(r, tq, noPivoting, epsilon) == std::min(m,n));
977 }
978 
979  /** Deprecated, use \ref linearSolveUpperTriangular().
980  */
981 template <class T, class C1, class C2, class C3>
982 inline
985 {
986  return linearSolveUpperTriangular(r, b, x);
987 }
988 
989  /** Solve a linear system with upper-triangular coefficient matrix.
990 
991  The square matrix \a r must be an upper-triangular coefficient matrix as can,
992  for example, be obtained by means of QR decomposition. If \a r doesn't have full rank
993  the function fails and returns <tt>false</tt>, otherwise it returns <tt>true</tt>. The
994  lower triangular part of matrix \a r will not be touched, so it doesn't need to contain zeros.
995 
996  The column vectors of matrix \a b are the right-hand sides of the equation (several equations
997  with the same coefficients can thus be solved in one go). The result is returned
998  int \a x, whose columns contain the solutions for the corresponding
999  columns of \a b. This implementation can be applied in-place, i.e. <tt>&b == &x</tt> is allowed.
1000  The following size requirements apply:
1001 
1002  \code
1003  rowCount(r) == columnCount(r);
1004  rowCount(r) == rowCount(b);
1005  columnCount(r) == rowCount(x);
1006  columnCount(b) == columnCount(x);
1007  \endcode
1008 
1009  <b>\#include</b> <vigra/linear_solve.hxx> or<br>
1010  <b>\#include</b> <vigra/linear_algebra.hxx><br>
1011  Namespaces: vigra and vigra::linalg
1012  */
1013 template <class T, class C1, class C2, class C3>
1016 {
1017  MultiArrayIndex m = rowCount(r);
1018  MultiArrayIndex rhsCount = columnCount(b);
1019  vigra_precondition(m == columnCount(r),
1020  "linearSolveUpperTriangular(): square coefficient matrix required.");
1021  vigra_precondition(m == rowCount(b) && m == rowCount(x) && rhsCount == columnCount(x),
1022  "linearSolveUpperTriangular(): matrix shape mismatch.");
1023 
1024  for(MultiArrayIndex k = 0; k < rhsCount; ++k)
1025  {
1026  for(int i=m-1; i>=0; --i)
1027  {
1028  if(r(i,i) == NumericTraits<T>::zero())
1029  return false; // r doesn't have full rank
1030  T sum = b(i, k);
1031  for(MultiArrayIndex j=i+1; j<m; ++j)
1032  sum -= r(i, j) * x(j, k);
1033  x(i, k) = sum / r(i, i);
1034  }
1035  }
1036  return true;
1037 }
1038 
1039  /** Solve a linear system with lower-triangular coefficient matrix.
1040 
1041  The square matrix \a l must be a lower-triangular coefficient matrix. If \a l
1042  doesn't have full rank the function fails and returns <tt>false</tt>,
1043  otherwise it returns <tt>true</tt>. The upper triangular part of matrix \a l will not be touched,
1044  so it doesn't need to contain zeros.
1045 
1046  The column vectors of matrix \a b are the right-hand sides of the equation (several equations
1047  with the same coefficients can thus be solved in one go). The result is returned
1048  in \a x, whose columns contain the solutions for the corresponding
1049  columns of \a b. This implementation can be applied in-place, i.e. <tt>&b == &x</tt> is allowed.
1050  The following size requirements apply:
1051 
1052  \code
1053  rowCount(l) == columnCount(l);
1054  rowCount(l) == rowCount(b);
1055  columnCount(l) == rowCount(x);
1056  columnCount(b) == columnCount(x);
1057  \endcode
1058 
1059  <b>\#include</b> <vigra/linear_solve.hxx> or<br>
1060  <b>\#include</b> <vigra/linear_algebra.hxx><br>
1061  Namespaces: vigra and vigra::linalg
1062  */
1063 template <class T, class C1, class C2, class C3>
1066 {
1069  vigra_precondition(m == rowCount(l),
1070  "linearSolveLowerTriangular(): square coefficient matrix required.");
1071  vigra_precondition(m == rowCount(b) && m == rowCount(x) && n == columnCount(x),
1072  "linearSolveLowerTriangular(): matrix shape mismatch.");
1073 
1074  for(MultiArrayIndex k = 0; k < n; ++k)
1075  {
1076  for(MultiArrayIndex i=0; i<m; ++i)
1077  {
1078  if(l(i,i) == NumericTraits<T>::zero())
1079  return false; // l doesn't have full rank
1080  T sum = b(i, k);
1081  for(MultiArrayIndex j=0; j<i; ++j)
1082  sum -= l(i, j) * x(j, k);
1083  x(i, k) = sum / l(i, i);
1084  }
1085  }
1086  return true;
1087 }
1088 
1089 
1090  /** Solve a linear system when the Cholesky decomposition of the left hand side is given.
1091 
1092  The square matrix \a L must be a lower-triangular matrix resulting from Cholesky
1093  decomposition of some positive definite coefficient matrix.
1094 
1095  The column vectors of matrix \a b are the right-hand sides of the equation (several equations
1096  with the same matrix \a L can thus be solved in one go). The result is returned
1097  in \a x, whose columns contain the solutions for the corresponding
1098  columns of \a b. This implementation can be applied in-place, i.e. <tt>&b == &x</tt> is allowed.
1099  The following size requirements apply:
1100 
1101  \code
1102  rowCount(L) == columnCount(L);
1103  rowCount(L) == rowCount(b);
1104  columnCount(L) == rowCount(x);
1105  columnCount(b) == columnCount(x);
1106  \endcode
1107 
1108  <b>\#include</b> <vigra/linear_solve.hxx> or<br>
1109  <b>\#include</b> <vigra/linear_algebra.hxx><br>
1110  Namespaces: vigra and vigra::linalg
1111  */
1112 template <class T, class C1, class C2, class C3>
1113 inline
1115 {
1116  /* Solve L * y = b */
1117  linearSolveLowerTriangular(L, b, x);
1118  /* Solve L^T * x = y */
1120 }
1121 
1122  /** Solve a linear system.
1123 
1124  <b> Declarations:</b>
1125 
1126  \code
1127  // use MultiArrayViews for input and output
1128  template <class T, class C1, class C2, class C3>
1129  bool linearSolve(MultiArrayView<2, T, C1> const & A,
1130  MultiArrayView<2, T, C2> const & b,
1131  MultiArrayView<2, T, C3> res,
1132  std::string method = "QR");
1133 
1134  // use TinyVector for RHS and result
1135  template <class T, class C1, int N>
1136  bool linearSolve(MultiArrayView<2, T, C1> const & A,
1137  TinyVector<T, N> const & b,
1138  TinyVector<T, N> & res,
1139  std::string method = "QR");
1140  \endcode
1141 
1142  \a A is the coefficient matrix, and the column vectors
1143  in \a b are the right-hand sides of the equation (so, several equations
1144  with the same coefficients can be solved in one go). The result is returned
1145  in \a res, whose columns contain the solutions for the corresponding
1146  columns of \a b. The number of columns of \a A must equal the number of rows of
1147  both \a b and \a res, and the number of columns of \a b and \a res must match.
1148  If right-hand-side and result are specified as TinyVector, the number of columns
1149  of \a A must equal N.
1150 
1151  \a method must be one of the following:
1152  <DL>
1153  <DT>"Cholesky"<DD> Compute the solution by means of Cholesky decomposition. The
1154  coefficient matrix \a A must by symmetric positive definite. If
1155  this is not the case, the function returns <tt>false</tt>.
1156 
1157  <DT>"QR"<DD> (default) Compute the solution by means of QR decomposition. The
1158  coefficient matrix \a A can be square or rectangular. In the latter case,
1159  it must have more rows than columns, and the solution will be computed in the
1160  least squares sense. If \a A doesn't have full rank, the function
1161  returns <tt>false</tt>.
1162 
1163  <DT>"SVD"<DD> Compute the solution by means of singular value decomposition. The
1164  coefficient matrix \a A can be square or rectangular. In the latter case,
1165  it must have more rows than columns, and the solution will be computed in the
1166  least squares sense. If \a A doesn't have full rank, the function
1167  returns <tt>false</tt>.
1168 
1169  <DT>"NE"<DD> Compute the solution by means of the normal equations, i.e. by applying Cholesky
1170  decomposition to the equivalent problem <tt>A'*A*x = A'*b</tt>. This only makes sense
1171  when the equation is to be solved in the least squares sense, i.e. when \a A is a
1172  rectangular matrix with more rows than columns. If \a A doesn't have full column rank,
1173  the function returns <tt>false</tt>.
1174  </DL>
1175 
1176  This function can be applied in-place, i.e. <tt>&b == &res</tt> or <tt>&A == &res</tt> are allowed
1177  (provided they have the required shapes).
1178 
1179  The following size requirements apply:
1180 
1181  \code
1182  rowCount(r) == rowCount(b);
1183  columnCount(r) == rowCount(x);
1184  columnCount(b) == columnCount(x);
1185  \endcode
1186 
1187  <b>\#include</b> <vigra/linear_solve.hxx> or<br>
1188  <b>\#include</b> <vigra/linear_algebra.hxx><br>
1189  Namespaces: vigra and vigra::linalg
1190  */
1191 doxygen_overloaded_function(template <...> bool linearSolve)
1192 
1193 template <class T, class C1, class C2, class C3>
1194 bool linearSolve(MultiArrayView<2, T, C1> const & A,
1195  MultiArrayView<2, T, C2> const & b,
1197  std::string method = "QR")
1198 {
1199  const MultiArrayIndex n = columnCount(A);
1200  const MultiArrayIndex m = rowCount(A);
1201 
1202  vigra_precondition(n <= m,
1203  "linearSolve(): Coefficient matrix A must have at least as many rows as columns.");
1204  vigra_precondition(n == rowCount(res) &&
1205  m == rowCount(b) && columnCount(b) == columnCount(res),
1206  "linearSolve(): matrix shape mismatch.");
1207 
1208  method = tolower(method);
1209  if(method == "cholesky")
1210  {
1211  vigra_precondition(columnCount(A) == rowCount(A),
1212  "linearSolve(): Cholesky method requires square coefficient matrix.");
1213  Matrix<T> L(A.shape());
1214  if(!choleskyDecomposition(A, L))
1215  return false; // false if A wasn't symmetric positive definite
1216  choleskySolve(L, b, res);
1217  }
1218  else if(method == "qr")
1219  {
1220  return (MultiArrayIndex)linearSolveQR(A, b, res) == n;
1221  }
1222  else if(method == "ne")
1223  {
1224  return linearSolve(transpose(A)*A, transpose(A)*b, res, "Cholesky");
1225  }
1226  else if(method == "svd")
1227  {
1228  MultiArrayIndex rhsCount = columnCount(b);
1229  Matrix<T> u(A.shape()), s(n, 1), v(n, n);
1230 
1232 
1233  Matrix<T> t = transpose(u)*b;
1234  for(MultiArrayIndex l=0; l<rhsCount; ++l)
1235  {
1236  for(MultiArrayIndex k=0; k<rank; ++k)
1237  t(k,l) /= s(k,0);
1238  for(MultiArrayIndex k=rank; k<n; ++k)
1239  t(k,l) = NumericTraits<T>::zero();
1240  }
1241  res = v*t;
1242 
1243  return (rank == n);
1244  }
1245  else
1246  {
1247  vigra_precondition(false, "linearSolve(): Unknown solution method.");
1248  }
1249  return true;
1250 }
1251 
1252 template <class T, class C1, int N>
1253 bool linearSolve(MultiArrayView<2, T, C1> const & A,
1254  TinyVector<T, N> const & b,
1255  TinyVector<T, N> & res,
1256  std::string method = "QR")
1257 {
1258  Shape2 shape(N, 1);
1259  return linearSolve(A, MultiArrayView<2, T>(shape, b.data()), MultiArrayView<2, T>(shape, res.data()), method);
1260 }
1261 
1262 //@}
1263 
1264 } // namespace linalg
1265 
1266 using linalg::inverse;
1267 using linalg::determinant;
1269 using linalg::linearSolve;
1270 using linalg::choleskySolve;
1275 
1276 } // namespace vigra
1277 
1278 
1279 #endif // VIGRA_LINEAR_SOLVE_HXX
bool qrDecomposition(MultiArrayView< 2, T, C1 > const &a, MultiArrayView< 2, T, C2 > &q, MultiArrayView< 2, T, C3 > &r, double epsilon=0.0)
Definition: linear_solve.hxx:962
FixedPoint16< 2, OverflowHandling > atan2(FixedPoint16< IntBits, OverflowHandling > y, FixedPoint16< IntBits, OverflowHandling > x)
Arctangent. Accuracy better than 1/3 degree (9 significant bits).
Definition: fixedpoint.hxx:1654
MultiArrayView< 2, T, C > columnVector(MultiArrayView< 2, T, C > const &m, MultiArrayIndex d)
Definition: matrix.hxx:725
MultiArrayIndex rowCount(const MultiArrayView< 2, T, C > &x)
Definition: matrix.hxx:669
void transpose(const MultiArrayView< 2, T, C1 > &v, MultiArrayView< 2, T, C2 > &r)
Definition: matrix.hxx:963
linalg::TemporaryMatrix< T > sin(MultiArrayView< 2, T, C > const &v)
bool linearSolveLowerTriangular(const MultiArrayView< 2, T, C1 > &l, const MultiArrayView< 2, T, C2 > &b, MultiArrayView< 2, T, C3 > x)
Definition: linear_solve.hxx:1064
const difference_type & shape() const
Definition: multi_array.hxx:1551
bool reverseElimination(const MultiArrayView< 2, T, C1 > &r, const MultiArrayView< 2, T, C2 > &b, MultiArrayView< 2, T, C3 > x)
Definition: linear_solve.hxx:983
bool isSymmetric(const MultiArrayView< 2, T, C > &v)
Definition: matrix.hxx:777
void choleskySolve(MultiArrayView< 2, T, C1 > const &L, MultiArrayView< 2, T, C2 > const &b, MultiArrayView< 2, T, C3 > &x)
Definition: linear_solve.hxx:1114
std::ptrdiff_t MultiArrayIndex
Definition: multi_shape.hxx:55
Definition: accessor.hxx:43
FixedPoint16< IntBits, OverflowHandling > hypot(FixedPoint16< IntBits, OverflowHandling > v1, FixedPoint16< IntBits, OverflowHandling > v2)
Length of hypotenuse.
Definition: fixedpoint.hxx:1640
NumericTraits< V >::Promote sum(TinyVectorBase< V, SIZE, D1, D2 > const &l)
sum of the vector's elements
Definition: tinyvector.hxx:1871
linalg::TemporaryMatrix< T > sq(MultiArrayView< 2, T, C > const &v)
Matrix< T, ALLLOC >::NormType norm(const Matrix< T, ALLLOC > &a)
Matrix< T, ALLLOC >::SquaredNormType squaredNorm(const Matrix< T, ALLLOC > &a)
Class for fixed size vectors.This class contains an array of size SIZE of the specified VALUETYPE...
Definition: accessor.hxx:940
int argMax(MultiArrayView< 2, T, C > const &a)
Find the index of the maximum element in a matrix.
Definition: matrix.hxx:1984
NormTraits< T >::SquaredNormType dot(const MultiArrayView< 2, T, C1 > &x, const MultiArrayView< 2, T, C2 > &y)
Definition: matrix.hxx:1340
double gamma(double x)
The gamma function.
Definition: mathutil.hxx:1549
MultiArrayShape< 2 >::type Shape2
shape type for MultiArray<2, T>
Definition: multi_shape.hxx:241
linalg::TemporaryMatrix< T > log(MultiArrayView< 2, T, C > const &v)
MultiArrayView< 2, T, C > rowVector(MultiArrayView< 2, T, C > const &m, MultiArrayIndex d)
Definition: matrix.hxx:695
MultiArrayIndex columnCount(const MultiArrayView< 2, T, C > &x)
Definition: matrix.hxx:682
bool choleskyDecomposition(MultiArrayView< 2, T, C1 > const &A, MultiArrayView< 2, T, C2 > &L)
Definition: linear_solve.hxx:907
T logDeterminant(MultiArrayView< 2, T, C1 > const &a)
Definition: linear_solve.hxx:859
bool linearSolveUpperTriangular(const MultiArrayView< 2, T, C1 > &r, const MultiArrayView< 2, T, C2 > &b, MultiArrayView< 2, T, C3 > x)
Definition: linear_solve.hxx:1014
unsigned int singularValueDecomposition(MultiArrayView< 2, T, C1 > const &A, MultiArrayView< 2, T, C2 > &U, MultiArrayView< 2, T, C3 > &S, MultiArrayView< 2, T, C4 > &V)
Definition: singular_value_decomposition.hxx:75
bool inverse(const MultiArrayView< 2, T, C1 > &v, MultiArrayView< 2, T, C2 > &res)
Definition: linear_solve.hxx:719
linalg::TemporaryMatrix< T > abs(MultiArrayView< 2, T, C > const &v)
linalg::TemporaryMatrix< T > cos(MultiArrayView< 2, T, C > const &v)
bool linearSolve(...)
SquareRootTraits< FixedPoint< IntBits, FracBits > >::SquareRootResult sqrt(FixedPoint< IntBits, FracBits > v)
square root.
Definition: fixedpoint.hxx:616
T determinant(MultiArrayView< 2, T, C1 > const &a, std::string method="LU")
Definition: linear_solve.hxx:816

© 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