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

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

© 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.11.1