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

eigensystem.hxx
1/************************************************************************/
2/* */
3/* Copyright 2004 by 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_EIGENSYSTEM_HXX
38#define VIGRA_EIGENSYSTEM_HXX
39
40#include <algorithm>
41#include <complex>
42#include "matrix.hxx"
43#include "array_vector.hxx"
44#include "polynomial.hxx"
45
46namespace vigra
47{
48
49namespace linalg
50{
51
52namespace detail
53{
54
55// code adapted from JAMA
56// a and b will be overwritten
57template <class T, class C1, class C2>
58void
59housholderTridiagonalization(MultiArrayView<2, T, C1> &a, MultiArrayView<2, T, C2> &b)
60{
61 const MultiArrayIndex n = rowCount(a);
62 vigra_precondition(n == columnCount(a),
63 "housholderTridiagonalization(): matrix must be square.");
64 vigra_precondition(n == rowCount(b) && 2 <= columnCount(b),
65 "housholderTridiagonalization(): matrix size mismatch.");
66
69
70 for(MultiArrayIndex j = 0; j < n; ++j)
71 {
72 d(j) = a(n-1, j);
73 }
74
75 // Householder reduction to tridiagonalMatrix form.
76
77 for(int i = n-1; i > 0; --i)
78 {
79 // Scale to avoid under/overflow.
80
81 T scale = 0.0;
82 T h = 0.0;
83 for(int k = 0; k < i; ++k)
84 {
85 scale = scale + abs(d(k));
86 }
87 if(scale == 0.0)
88 {
89 e(i) = d(i-1);
90 for(int j = 0; j < i; ++j)
91 {
92 d(j) = a(i-1, j);
93 a(i, j) = 0.0;
94 a(j, i) = 0.0;
95 }
96 }
97 else
98 {
99 // Generate Householder vector.
100
101 for(int k = 0; k < i; ++k)
102 {
103 d(k) /= scale;
104 h += sq(d(k));
105 }
106 T f = d(i-1);
107 T g = VIGRA_CSTD::sqrt(h);
108 if(f > 0) {
109 g = -g;
110 }
111 e(i) = scale * g;
112 h -= f * g;
113 d(i-1) = f - g;
114 for(int j = 0; j < i; ++j)
115 {
116 e(j) = 0.0;
117 }
118
119 // Apply similarity transformation to remaining columns.
120
121 for(int j = 0; j < i; ++j)
122 {
123 f = d(j);
124 a(j, i) = f;
125 g = e(j) + a(j, j) * f;
126 for(int k = j+1; k <= i-1; ++k)
127 {
128 g += a(k, j) * d(k);
129 e(k) += a(k, j) * f;
130 }
131 e(j) = g;
132 }
133 f = 0.0;
134 for(int j = 0; j < i; ++j)
135 {
136 e(j) /= h;
137 f += e(j) * d(j);
138 }
139 T hh = f / (h + h);
140 for(int j = 0; j < i; ++j)
141 {
142 e(j) -= hh * d(j);
143 }
144 for(int j = 0; j < i; ++j)
145 {
146 f = d(j);
147 g = e(j);
148 for(int k = j; k <= i-1; ++k)
149 {
150 a(k, j) -= (f * e(k) + g * d(k));
151 }
152 d(j) = a(i-1, j);
153 a(i, j) = 0.0;
154 }
155 }
156 d(i) = h;
157 }
158
159 // Accumulate transformations.
160
161 for(MultiArrayIndex i = 0; i < n-1; ++i)
162 {
163 a(n-1, i) = a(i, i);
164 a(i, i) = 1.0;
165 T h = d(i+1);
166 if(h != 0.0)
167 {
168 for(MultiArrayIndex k = 0; k <= i; ++k)
169 {
170 d(k) = a(k, i+1) / h;
171 }
172 for(MultiArrayIndex j = 0; j <= i; ++j)
173 {
174 T g = 0.0;
175 for(MultiArrayIndex k = 0; k <= i; ++k)
176 {
177 g += a(k, i+1) * a(k, j);
178 }
179 for(MultiArrayIndex k = 0; k <= i; ++k)
180 {
181 a(k, j) -= g * d(k);
182 }
183 }
184 }
185 for(MultiArrayIndex k = 0; k <= i; ++k)
186 {
187 a(k, i+1) = 0.0;
188 }
189 }
190 for(MultiArrayIndex j = 0; j < n; ++j)
191 {
192 d(j) = a(n-1, j);
193 a(n-1, j) = 0.0;
194 }
195 a(n-1, n-1) = 1.0;
196 e(0) = 0.0;
197}
198
199// code adapted from JAMA
200// de and z will be overwritten
201template <class T, class C1, class C2>
202bool
203tridiagonalMatrixEigensystem(MultiArrayView<2, T, C1> &de, MultiArrayView<2, T, C2> &z)
204{
206 vigra_precondition(n == columnCount(z),
207 "tridiagonalMatrixEigensystem(): matrix must be square.");
208 vigra_precondition(n == rowCount(de) && 2 <= columnCount(de),
209 "tridiagonalMatrixEigensystem(): matrix size mismatch.");
210
213
214 for(MultiArrayIndex i = 1; i < n; i++) {
215 e(i-1) = e(i);
216 }
217 e(n-1) = 0.0;
218
219 T f = 0.0;
220 T tst1 = 0.0;
221 T eps = VIGRA_CSTD::pow(2.0,-52.0);
222 for(MultiArrayIndex l = 0; l < n; ++l)
223 {
224 // Find small subdiagonalMatrix element
225
226 tst1 = std::max(tst1, abs(d(l)) + abs(e(l)));
227 MultiArrayIndex m = l;
228
229 // Original while-loop from Java code
230 while(m < n)
231 {
232 if(abs(e(m)) <= eps*tst1)
233 break;
234 ++m;
235 }
236
237 // If m == l, d(l) is an eigenvalue,
238 // otherwise, iterate.
239
240 if(m > l)
241 {
242 int iter = 0;
243 do
244 {
245 if(++iter > 50)
246 return false; // too many iterations
247
248 // Compute implicit shift
249
250 T g = d(l);
251 T p = (d(l+1) - g) / (2.0 * e(l));
252 T r = hypot(p,1.0);
253 if(p < 0)
254 {
255 r = -r;
256 }
257 d(l) = e(l) / (p + r);
258 d(l+1) = e(l) * (p + r);
259 T dl1 = d(l+1);
260 T h = g - d(l);
261 for(MultiArrayIndex i = l+2; i < n; ++i)
262 {
263 d(i) -= h;
264 }
265 f = f + h;
266
267 // Implicit QL transformation.
268
269 p = d(m);
270 T c = 1.0;
271 T c2 = c;
272 T c3 = c;
273 T el1 = e(l+1);
274 T s = 0.0;
275 T s2 = 0.0;
276 for(int i = m-1; i >= (int)l; --i)
277 {
278 c3 = c2;
279 c2 = c;
280 s2 = s;
281 g = c * e(i);
282 h = c * p;
283 r = hypot(p,e(i));
284 e(i+1) = s * r;
285 s = e(i) / r;
286 c = p / r;
287 p = c * d(i) - s * g;
288 d(i+1) = h + s * (c * g + s * d(i));
289
290 // Accumulate transformation.
291
292 for(MultiArrayIndex k = 0; k < n; ++k)
293 {
294 h = z(k, i+1);
295 z(k, i+1) = s * z(k, i) + c * h;
296 z(k, i) = c * z(k, i) - s * h;
297 }
298 }
299 p = -s * s2 * c3 * el1 * e(l) / dl1;
300 e(l) = s * p;
301 d(l) = c * p;
302
303 // Check for convergence.
304
305 } while(abs(e(l)) > eps*tst1);
306 }
307 d(l) = d(l) + f;
308 e(l) = 0.0;
309 }
310
311 // Sort eigenvalues and corresponding vectors.
312
313 for(MultiArrayIndex i = 0; i < n-1; ++i)
314 {
315 MultiArrayIndex k = i;
316 T p = d(i);
317 for(MultiArrayIndex j = i+1; j < n; ++j)
318 {
319 T p1 = d(j);
320 if(p < p1)
321 {
322 k = j;
323 p = p1;
324 }
325 }
326 if(k != i)
327 {
328 std::swap(d(k), d(i));
329 for(MultiArrayIndex j = 0; j < n; ++j)
330 {
331 std::swap(z(j, i), z(j, k));
332 }
333 }
334 }
335 return true;
336}
337
338// Nonsymmetric reduction to Hessenberg form.
339
340template <class T, class C1, class C2>
341void nonsymmetricHessenbergReduction(MultiArrayView<2, T, C1> & H, MultiArrayView<2, T, C2> & V)
342{
343 // This is derived from the Algol procedures orthes and ortran,
344 // by Martin and Wilkinson, Handbook for Auto. Comp.,
345 // Vol.ii-Linear Algebra, and the corresponding
346 // Fortran subroutines in EISPACK.
347
348 int n = rowCount(H);
349 int low = 0;
350 int high = n-1;
351 ArrayVector<T> ort(n);
352
353 for(int m = low+1; m <= high-1; ++m)
354 {
355 // Scale column.
356
357 T scale = 0.0;
358 for(int i = m; i <= high; ++i)
359 {
360 scale = scale + abs(H(i, m-1));
361 }
362 if(scale != 0.0)
363 {
364
365 // Compute Householder transformation.
366
367 T h = 0.0;
368 for(int i = high; i >= m; --i)
369 {
370 ort[i] = H(i, m-1)/scale;
371 h += sq(ort[i]);
372 }
373 T g = VIGRA_CSTD::sqrt(h);
374 if(ort[m] > 0)
375 {
376 g = -g;
377 }
378 h = h - ort[m] * g;
379 ort[m] = ort[m] - g;
380
381 // Apply Householder similarity transformation
382 // H = (I-u*u'/h)*H*(I-u*u')/h)
383
384 for(int j = m; j < n; ++j)
385 {
386 T f = 0.0;
387 for(int i = high; i >= m; --i)
388 {
389 f += ort[i]*H(i, j);
390 }
391 f = f/h;
392 for(int i = m; i <= high; ++i)
393 {
394 H(i, j) -= f*ort[i];
395 }
396 }
397
398 for(int i = 0; i <= high; ++i)
399 {
400 T f = 0.0;
401 for(int j = high; j >= m; --j)
402 {
403 f += ort[j]*H(i, j);
404 }
405 f = f/h;
406 for(int j = m; j <= high; ++j)
407 {
408 H(i, j) -= f*ort[j];
409 }
410 }
411 ort[m] = scale*ort[m];
412 H(m, m-1) = scale*g;
413 }
414 }
415
416 // Accumulate transformations (Algol's ortran).
417
418 for(int i = 0; i < n; ++i)
419 {
420 for(int j = 0; j < n; ++j)
421 {
422 V(i, j) = (i == j ? 1.0 : 0.0);
423 }
424 }
425
426 for(int m = high-1; m >= low+1; --m)
427 {
428 if(H(m, m-1) != 0.0)
429 {
430 for(int i = m+1; i <= high; ++i)
431 {
432 ort[i] = H(i, m-1);
433 }
434 for(int j = m; j <= high; ++j)
435 {
436 T g = 0.0;
437 for(int i = m; i <= high; ++i)
438 {
439 g += ort[i] * V(i, j);
440 }
441 // Double division avoids possible underflow
442 g = (g / ort[m]) / H(m, m-1);
443 for(int i = m; i <= high; ++i)
444 {
445 V(i, j) += g * ort[i];
446 }
447 }
448 }
449 }
450}
451
452
453// Complex scalar division.
454
455template <class T>
456void cdiv(T xr, T xi, T yr, T yi, T & cdivr, T & cdivi)
457{
458 T r,d;
459 if(abs(yr) > abs(yi))
460 {
461 r = yi/yr;
462 d = yr + r*yi;
463 cdivr = (xr + r*xi)/d;
464 cdivi = (xi - r*xr)/d;
465 }
466 else
467 {
468 r = yr/yi;
469 d = yi + r*yr;
470 cdivr = (r*xr + xi)/d;
471 cdivi = (r*xi - xr)/d;
472 }
473}
474
475template <class T, class C>
476int hessenbergQrDecompositionHelper(MultiArrayView<2, T, C> & H, int n, int l, double eps,
477 T & p, T & q, T & r, T & s, T & w, T & x, T & y, T & z)
478{
479 int m = n-2;
480 while(m >= l)
481 {
482 z = H(m, m);
483 r = x - z;
484 s = y - z;
485 p = (r * s - w) / H(m+1, m) + H(m, m+1);
486 q = H(m+1, m+1) - z - r - s;
487 r = H(m+2, m+1);
488 s = abs(p) + abs(q) + abs(r);
489 p = p / s;
490 q = q / s;
491 r = r / s;
492 if(m == l)
493 {
494 break;
495 }
496 if(abs(H(m, m-1)) * (abs(q) + abs(r)) <
497 eps * (abs(p) * (abs(H(m-1, m-1)) + abs(z) +
498 abs(H(m+1, m+1)))))
499 {
500 break;
501 }
502 --m;
503 }
504 return m;
505}
506
507
508
509// Nonsymmetric reduction from Hessenberg to real Schur form.
510
511template <class T, class C1, class C2, class C3>
512bool hessenbergQrDecomposition(MultiArrayView<2, T, C1> & H, MultiArrayView<2, T, C2> & V,
514{
515
516 // This is derived from the Algol procedure hqr2,
517 // by Martin and Wilkinson, Handbook for Auto. Comp.,
518 // Vol.ii-Linear Algebra, and the corresponding
519 // Fortran subroutine in EISPACK.
520
521 // Initialize
524
525 int nn = rowCount(H);
526 int n = nn-1;
527 int low = 0;
528 int high = nn-1;
529 T eps = VIGRA_CSTD::pow(2.0, sizeof(T) == sizeof(float)
530 ? -23.0
531 : -52.0);
532 T exshift = 0.0;
533 T p=0,q=0,r=0,s=0,z=0,t,w,x,y;
534 T norm = vigra::norm(H);
535
536 // Outer loop over eigenvalue index
537 int iter = 0;
538 while(n >= low)
539 {
540
541 // Look for single small sub-diagonal element
542 int l = n;
543 while (l > low)
544 {
545 s = abs(H(l-1, l-1)) + abs(H(l, l));
546 if(s == 0.0)
547 {
548 s = norm;
549 }
550 if(abs(H(l, l-1)) < eps * s)
551 {
552 break;
553 }
554 --l;
555 }
556
557 // Check for convergence
558 // One root found
559 if(l == n)
560 {
561 H(n, n) = H(n, n) + exshift;
562 d(n) = H(n, n);
563 e(n) = 0.0;
564 --n;
565 iter = 0;
566
567 // Two roots found
568
569 }
570 else if(l == n-1)
571 {
572 w = H(n, n-1) * H(n-1, n);
573 p = (H(n-1, n-1) - H(n, n)) / 2.0;
574 q = p * p + w;
575 z = VIGRA_CSTD::sqrt(abs(q));
576 H(n, n) = H(n, n) + exshift;
577 H(n-1, n-1) = H(n-1, n-1) + exshift;
578 x = H(n, n);
579
580 // Real pair
581
582 if(q >= 0)
583 {
584 if(p >= 0)
585 {
586 z = p + z;
587 }
588 else
589 {
590 z = p - z;
591 }
592 d(n-1) = x + z;
593 d(n) = d(n-1);
594 if(z != 0.0)
595 {
596 d(n) = x - w / z;
597 }
598 e(n-1) = 0.0;
599 e(n) = 0.0;
600 x = H(n, n-1);
601 s = abs(x) + abs(z);
602 p = x / s;
603 q = z / s;
604 r = VIGRA_CSTD::sqrt(p * p+q * q);
605 p = p / r;
606 q = q / r;
607
608 // Row modification
609
610 for(int j = n-1; j < nn; ++j)
611 {
612 z = H(n-1, j);
613 H(n-1, j) = q * z + p * H(n, j);
614 H(n, j) = q * H(n, j) - p * z;
615 }
616
617 // Column modification
618
619 for(int i = 0; i <= n; ++i)
620 {
621 z = H(i, n-1);
622 H(i, n-1) = q * z + p * H(i, n);
623 H(i, n) = q * H(i, n) - p * z;
624 }
625
626 // Accumulate transformations
627
628 for(int i = low; i <= high; ++i)
629 {
630 z = V(i, n-1);
631 V(i, n-1) = q * z + p * V(i, n);
632 V(i, n) = q * V(i, n) - p * z;
633 }
634
635 // Complex pair
636
637 }
638 else
639 {
640 d(n-1) = x + p;
641 d(n) = x + p;
642 e(n-1) = z;
643 e(n) = -z;
644 }
645 n = n - 2;
646 iter = 0;
647
648 // No convergence yet
649
650 }
651 else
652 {
653
654 // Form shift
655
656 x = H(n, n);
657 y = 0.0;
658 w = 0.0;
659 if(l < n)
660 {
661 y = H(n-1, n-1);
662 w = H(n, n-1) * H(n-1, n);
663 }
664
665 // Wilkinson's original ad hoc shift
666
667 if(iter == 10)
668 {
669 exshift += x;
670 for(int i = low; i <= n; ++i)
671 {
672 H(i, i) -= x;
673 }
674 s = abs(H(n, n-1)) + abs(H(n-1, n-2));
675 x = y = 0.75 * s;
676 w = -0.4375 * s * s;
677 }
678
679 // MATLAB's new ad hoc shift
680
681 if(iter == 30)
682 {
683 s = (y - x) / 2.0;
684 s = s * s + w;
685 if(s > 0)
686 {
687 s = VIGRA_CSTD::sqrt(s);
688 if(y < x)
689 {
690 s = -s;
691 }
692 s = x - w / ((y - x) / 2.0 + s);
693 for(int i = low; i <= n; ++i)
694 {
695 H(i, i) -= s;
696 }
697 exshift += s;
698 x = y = w = 0.964;
699 }
700 }
701
702 iter = iter + 1;
703 if(iter > 60)
704 return false;
705
706 // Look for two consecutive small sub-diagonal elements
707 int m = hessenbergQrDecompositionHelper(H, n, l, eps, p, q, r, s, w, x, y, z);
708 for(int i = m+2; i <= n; ++i)
709 {
710 H(i, i-2) = 0.0;
711 if(i > m+2)
712 {
713 H(i, i-3) = 0.0;
714 }
715 }
716
717 // Double QR step involving rows l:n and columns m:n
718
719 for(int k = m; k <= n-1; ++k)
720 {
721 int notlast = (k != n-1);
722 if(k != m) {
723 p = H(k, k-1);
724 q = H(k+1, k-1);
725 r = (notlast ? H(k+2, k-1) : 0.0);
726 x = abs(p) + abs(q) + abs(r);
727 if(x != 0.0)
728 {
729 p = p / x;
730 q = q / x;
731 r = r / x;
732 }
733 }
734 if(x == 0.0)
735 {
736 break;
737 }
738 s = VIGRA_CSTD::sqrt(p * p + q * q + r * r);
739 if(p < 0)
740 {
741 s = -s;
742 }
743 if(s != 0)
744 {
745 if(k != m)
746 {
747 H(k, k-1) = -s * x;
748 }
749 else if(l != m)
750 {
751 H(k, k-1) = -H(k, k-1);
752 }
753 p = p + s;
754 x = p / s;
755 y = q / s;
756 z = r / s;
757 q = q / p;
758 r = r / p;
759
760 // Row modification
761
762 for(int j = k; j < nn; ++j)
763 {
764 p = H(k, j) + q * H(k+1, j);
765 if(notlast)
766 {
767 p = p + r * H(k+2, j);
768 H(k+2, j) = H(k+2, j) - p * z;
769 }
770 H(k, j) = H(k, j) - p * x;
771 H(k+1, j) = H(k+1, j) - p * y;
772 }
773
774 // Column modification
775
776 for(int i = 0; i <= std::min(n,k+3); ++i)
777 {
778 p = x * H(i, k) + y * H(i, k+1);
779 if(notlast)
780 {
781 p = p + z * H(i, k+2);
782 H(i, k+2) = H(i, k+2) - p * r;
783 }
784 H(i, k) = H(i, k) - p;
785 H(i, k+1) = H(i, k+1) - p * q;
786 }
787
788 // Accumulate transformations
789
790 for(int i = low; i <= high; ++i)
791 {
792 p = x * V(i, k) + y * V(i, k+1);
793 if(notlast)
794 {
795 p = p + z * V(i, k+2);
796 V(i, k+2) = V(i, k+2) - p * r;
797 }
798 V(i, k) = V(i, k) - p;
799 V(i, k+1) = V(i, k+1) - p * q;
800 }
801 } // (s != 0)
802 } // k loop
803 } // check convergence
804 } // while (n >= low)
805
806 // Backsubstitute to find vectors of upper triangular form
807
808 if(norm == 0.0)
809 {
810 return false;
811 }
812
813 for(n = nn-1; n >= 0; --n)
814 {
815 p = d(n);
816 q = e(n);
817
818 // Real vector
819
820 if(q == 0)
821 {
822 int l = n;
823 H(n, n) = 1.0;
824 for(int i = n-1; i >= 0; --i)
825 {
826 w = H(i, i) - p;
827 r = 0.0;
828 for(int j = l; j <= n; ++j)
829 {
830 r = r + H(i, j) * H(j, n);
831 }
832 if(e(i) < 0.0)
833 {
834 z = w;
835 s = r;
836 }
837 else
838 {
839 l = i;
840 if(e(i) == 0.0)
841 {
842 if(w != 0.0)
843 {
844 H(i, n) = -r / w;
845 }
846 else
847 {
848 H(i, n) = -r / (eps * norm);
849 }
850
851 // Solve real equations
852
853 }
854 else
855 {
856 x = H(i, i+1);
857 y = H(i+1, i);
858 q = (d(i) - p) * (d(i) - p) + e(i) * e(i);
859 t = (x * s - z * r) / q;
860 H(i, n) = t;
861 if(abs(x) > abs(z))
862 {
863 H(i+1, n) = (-r - w * t) / x;
864 }
865 else
866 {
867 H(i+1, n) = (-s - y * t) / z;
868 }
869 }
870
871 // Overflow control
872
873 t = abs(H(i, n));
874 if((eps * t) * t > 1)
875 {
876 for(int j = i; j <= n; ++j)
877 {
878 H(j, n) = H(j, n) / t;
879 }
880 }
881 }
882 }
883
884 // Complex vector
885
886 }
887 else if(q < 0)
888 {
889 int l = n-1;
890
891 // Last vector component imaginary so matrix is triangular
892
893 if(abs(H(n, n-1)) > abs(H(n-1, n)))
894 {
895 H(n-1, n-1) = q / H(n, n-1);
896 H(n-1, n) = -(H(n, n) - p) / H(n, n-1);
897 }
898 else
899 {
900 cdiv(0.0,-H(n-1, n),H(n-1, n-1)-p,q, H(n-1, n-1), H(n-1, n));
901 }
902 H(n, n-1) = 0.0;
903 H(n, n) = 1.0;
904 for(int i = n-2; i >= 0; --i)
905 {
906 T ra,sa,vr,vi;
907 ra = 0.0;
908 sa = 0.0;
909 for(int j = l; j <= n; ++j)
910 {
911 ra = ra + H(j, n-1) * H(i, j);
912 sa = sa + H(j, n) * H(i, j);
913 }
914 w = H(i, i) - p;
915
916 if(e(i) < 0.0)
917 {
918 z = w;
919 r = ra;
920 s = sa;
921 }
922 else
923 {
924 l = i;
925 if(e(i) == 0)
926 {
927 cdiv(-ra,-sa,w,q, H(i, n-1), H(i, n));
928 }
929 else
930 {
931 // Solve complex equations
932
933 x = H(i, i+1);
934 y = H(i+1, i);
935 vr = (d(i) - p) * (d(i) - p) + e(i) * e(i) - q * q;
936 vi = (d(i) - p) * 2.0 * q;
937 if((vr == 0.0) && (vi == 0.0))
938 {
939 vr = eps * norm * (abs(w) + abs(q) +
940 abs(x) + abs(y) + abs(z));
941 }
942 cdiv(x*r-z*ra+q*sa,x*s-z*sa-q*ra,vr,vi, H(i, n-1), H(i, n));
943 if(abs(x) > (abs(z) + abs(q)))
944 {
945 H(i+1, n-1) = (-ra - w * H(i, n-1) + q * H(i, n)) / x;
946 H(i+1, n) = (-sa - w * H(i, n) - q * H(i, n-1)) / x;
947 }
948 else
949 {
950 cdiv(-r-y*H(i, n-1),-s-y*H(i, n),z,q, H(i+1, n-1), H(i+1, n));
951 }
952 }
953
954 // Overflow control
955
956 t = std::max(abs(H(i, n-1)),abs(H(i, n)));
957 if((eps * t) * t > 1)
958 {
959 for(int j = i; j <= n; ++j)
960 {
961 H(j, n-1) = H(j, n-1) / t;
962 H(j, n) = H(j, n) / t;
963 }
964 }
965 }
966 }
967 }
968 }
969
970 // Back transformation to get eigenvectors of original matrix
971
972 for(int j = nn-1; j >= low; --j)
973 {
974 for(int i = low; i <= high; ++i)
975 {
976 z = 0.0;
977 for(int k = low; k <= std::min(j,high); ++k)
978 {
979 z = z + V(i, k) * H(k, j);
980 }
981 V(i, j) = z;
982 }
983 }
984 return true;
985}
986
987} // namespace detail
988
989/** \addtogroup MatrixAlgebra
990*/
991//@{
992 /** Compute the eigensystem of a symmetric matrix.
993
994 \a a is a real symmetric matrix, \a ew is a single-column matrix
995 holding the eigenvalues, and \a ev is a matrix of the same size as
996 \a a whose columns are the corresponding eigenvectors. Eigenvalues
997 will be sorted from largest to smallest.
998 The algorithm returns <tt>false</tt> when it doesn't
999 converge. It can be applied in-place, i.e. <tt>&a == &ev</tt> is allowed.
1000 The code of this function was adapted from JAMA.
1001
1002 <b>\#include</b> <vigra/eigensystem.hxx> or<br>
1003 <b>\#include</b> <vigra/linear_algebra.hxx><br>
1004 Namespaces: vigra and vigra::linalg
1005 */
1006template <class T, class C1, class C2, class C3>
1007bool
1010{
1011 vigra_precondition(isSymmetric(a),
1012 "symmetricEigensystem(): symmetric input matrix required.");
1013 const MultiArrayIndex acols = columnCount(a);
1014 vigra_precondition(1 == columnCount(ew) && acols == rowCount(ew) &&
1015 acols == columnCount(ev) && acols == rowCount(ev),
1016 "symmetricEigensystem(): matrix shape mismatch.");
1017
1018 ev.copy(a); // does nothing if &ev == &a
1019 Matrix<T> de(acols, 2);
1020 detail::housholderTridiagonalization(ev, de);
1021 if(!detail::tridiagonalMatrixEigensystem(de, ev))
1022 return false;
1023
1024 ew.copy(columnVector(de, 0));
1025 return true;
1026}
1027
1028namespace detail{
1029
1030template <class T, class C2, class C3>
1031bool
1032symmetricEigensystem2x2(T a00, T a01, T a11,
1034{
1035 double evec[2]={0,0};
1036
1037 /* Eigenvectors*/
1038 if (a01==0){
1039 if (fabs(a11)>fabs(a00)){
1040 evec[0]=0.;
1041 evec[1]=1.;
1042 ew(0,0)=a11;
1043 ew(1,0)=a00;
1044 }
1045 else if(fabs(a00)>fabs(a11)) {
1046 evec[0]=1.;
1047 evec[1]=0.;
1048 ew(0,0)=a00;
1049 ew(1,0)=a11;
1050 }
1051 else {
1052 evec[0]=.5* M_SQRT2;
1053 evec[1]=.5* M_SQRT2;
1054 ew(0,0)=a00;
1055 ew(1,0)=a11;
1056 }
1057 }
1058 else{
1059 double temp=a11-a00;
1060
1061 double coherence=sqrt(temp*temp+4*a01*a01);
1062 evec[0]=2*a01;
1063 evec[1]=temp+coherence;
1064 temp=std::sqrt(evec[0]*evec[0]+evec[1]*evec[1]);
1065 if (temp==0){
1066 evec[0]=.5* M_SQRT2;
1067 evec[1]=.5* M_SQRT2;
1068 ew(0,0)=1.;
1069 ew(1,0)=1.;
1070 }
1071 else{
1072 evec[0]/=temp;
1073 evec[1]/=temp;
1074
1075 /* Eigenvalues */
1076 ew(0,0)=.5*(a00+a11+coherence);
1077 ew(1,0)=.5*(a00+a11-coherence);
1078 }
1079 }
1080 ev(0,0)= evec[0];
1081 ev(1,0)= evec[1];
1082 ev(0,1)=-evec[1];
1083 ev(1,1)= evec[0];
1084 return true;
1085}
1086
1087template <class T, class C2, class C3>
1088bool
1089symmetricEigensystem3x3(T a00, T a01, T a02, T a11, T a12, T a22,
1090 MultiArrayView<2, T, C2> & ew, MultiArrayView<2, T, C3> & ev)
1091{
1092 symmetric3x3Eigenvalues(a00, a01, a02, a11, a12, a22,
1093 &ew(0,0), &ew(1,0), &ew(2,0));
1094
1095 /* Calculate eigen vectors */
1096 double a1=a01*a12,
1097 a2=a01*a02,
1098 a3=sq(a01);
1099
1100 double b1=a00-ew(0,0),
1101 b2=a11-ew(0,0);
1102 ev(0,0)=a1-a02*b2;
1103 ev(1,0)=a2-a12*b1;
1104 ev(2,0)=b1*b2-a3;
1105
1106 b1=a00-ew(1,0);
1107 b2=a11-ew(1,0);
1108 ev(0,1)=a1-a02*b2;
1109 ev(1,1)=a2-a12*b1;
1110 ev(2,1)=b1*b2-a3;
1111
1112 b1=a00-ew(2,0);
1113 b2=a11-ew(2,0);
1114 ev(0,2)=a1-a02*b2;
1115 ev(1,2)=a2-a12*b1;
1116 ev(2,2)=b1*b2-a3;
1117
1118 /* Eigen vector normalization */
1119 double l0=norm(columnVector(ev, 0));
1120 double l1=norm(columnVector(ev, 1));
1121 double l2=norm(columnVector(ev, 2));
1122
1123 /* Detect fail : eigenvectors with only zeros */
1124 double M = std::max(std::max(abs(ew(0,0)), abs(ew(1,0))), abs(ew(2,0)));
1125 double epsilon = 1e-12*M;
1126 if(l0<epsilon) { return false; }
1127 if(l1<epsilon) { return false; }
1128 if(l2<epsilon) { return false; }
1129
1130 columnVector(ev, 0) /= l0;
1131 columnVector(ev, 1) /= l1;
1132 columnVector(ev, 2) /= l2;
1133
1134 /* Succes */
1135 return true;
1136}
1137
1138} // closing namespace detail
1139
1140
1141 /** Fast computation of the eigensystem of a 2x2 or 3x3 symmetric matrix.
1142
1143 The function works like \ref symmetricEigensystem(), but uses fast analytic
1144 formula to avoid iterative computations.
1145
1146 <b>\#include</b> <vigra/eigensystem.hxx> or<br>
1147 <b>\#include</b> <vigra/linear_algebra.hxx><br>
1148 Namespaces: vigra and vigra::linalg
1149 */
1150template <class T, class C1, class C2, class C3>
1151bool
1154{
1155 vigra_precondition(isSymmetric(a),
1156 "symmetricEigensystemNoniterative(): symmetric input matrix required.");
1157 const MultiArrayIndex acols = columnCount(a);
1158 if(acols == 2)
1159 {
1160 detail::symmetricEigensystem2x2(a(0,0), a(0,1), a(1,1), ew, ev);
1161 return true;
1162 }
1163 if(acols == 3)
1164 {
1165 // try the fast algorithm
1166 if(detail::symmetricEigensystem3x3(a(0,0), a(0,1), a(0,2), a(1,1), a(1,2), a(2,2),
1167 ew, ev))
1168 return true;
1169 // fast algorithm failed => fall-back to iterative algorithm
1170 return symmetricEigensystem(a, ew, ev);
1171 }
1172 vigra_precondition(false,
1173 "symmetricEigensystemNoniterative(): can only handle 2x2 and 3x3 matrices.");
1174 return false;
1175}
1176
1177 /** Compute the eigensystem of a square, but
1178 not necessarily symmetric matrix.
1179
1180 \a a is a real square matrix, \a ew is a single-column matrix
1181 holding the possibly complex eigenvalues, and \a ev is a matrix of
1182 the same size as \a a whose columns are the corresponding eigenvectors.
1183 Eigenvalues will be sorted from largest to smallest magnitude.
1184 The algorithm returns <tt>false</tt> when it doesn't
1185 converge. It can be applied in-place, i.e. <tt>&a == &ev</tt> is allowed.
1186 The code of this function was adapted from JAMA.
1187
1188 <b>\#include</b> <vigra/eigensystem.hxx> or<br>
1189 <b>\#include</b> <vigra/linear_algebra.hxx><br>
1190 Namespaces: vigra and vigra::linalg
1191 */
1192template <class T, class C1, class C2, class C3>
1193bool
1195 MultiArrayView<2, std::complex<T>, C2> & ew, MultiArrayView<2, T, C3> & ev)
1196{
1197 const MultiArrayIndex acols = columnCount(a);
1198 vigra_precondition(acols == rowCount(a),
1199 "nonsymmetricEigensystem(): square input matrix required.");
1200 vigra_precondition(1 == columnCount(ew) && acols == rowCount(ew) &&
1201 acols == columnCount(ev) && acols == rowCount(ev),
1202 "nonsymmetricEigensystem(): matrix shape mismatch.");
1203
1204 Matrix<T> H(a);
1205 Matrix<T> de(acols, 2);
1206 detail::nonsymmetricHessenbergReduction(H, ev);
1207 if(!detail::hessenbergQrDecomposition(H, ev, de))
1208 return false;
1209
1210 for(MultiArrayIndex i = 0; i < acols; ++i)
1211 {
1212 ew(i,0) = std::complex<T>(de(i, 0), de(i, 1));
1213 }
1214 return true;
1215}
1216
1217 /** Compute the roots of a polynomial using the eigenvalue method.
1218
1219 \a poly is a real polynomial (compatible to \ref vigra::PolynomialView),
1220 and \a roots a complex valued vector (compatible to <tt>std::vector</tt>
1221 with a <tt>value_type</tt> compatible to the type <tt>POLYNOMIAL::Complex</tt>) to which
1222 the roots are appended. The function calls \ref nonsymmetricEigensystem() with the standard
1223 companion matrix yielding the roots as eigenvalues. It returns <tt>false</tt> if
1224 it fails to converge.
1225
1226 <b>\#include</b> <vigra/eigensystem.hxx> or<br>
1227 <b>\#include</b> <vigra/linear_algebra.hxx><br>
1228 Namespaces: vigra and vigra::linalg
1229
1230 \see polynomialRoots(), vigra::Polynomial
1231 */
1232template <class POLYNOMIAL, class VECTOR>
1233bool polynomialRootsEigenvalueMethod(POLYNOMIAL const & poly, VECTOR & roots, bool polishRoots)
1234{
1235 typedef typename POLYNOMIAL::value_type T;
1236 typedef typename POLYNOMIAL::Real Real;
1237 typedef typename POLYNOMIAL::Complex Complex;
1238 typedef Matrix<T> TMatrix;
1239 typedef Matrix<Complex> ComplexMatrix;
1240
1241 int const degree = poly.order();
1242 double const eps = poly.epsilon();
1243
1244 TMatrix inMatrix(degree, degree);
1245 for(int i = 0; i < degree; ++i)
1246 inMatrix(0, i) = -poly[degree - i - 1] / poly[degree];
1247 for(int i = 0; i < degree - 1; ++i)
1248 inMatrix(i + 1, i) = NumericTraits<T>::one();
1249 ComplexMatrix ew(degree, 1);
1250 TMatrix ev(degree, degree);
1251 bool success = nonsymmetricEigensystem(inMatrix, ew, ev);
1252 if(!success)
1253 return false;
1254 for(int i = 0; i < degree; ++i)
1255 {
1256 if(polishRoots)
1257 vigra::detail::laguerre1Root(poly, ew(i,0), 1);
1258 roots.push_back(vigra::detail::deleteBelowEpsilon(ew(i,0), eps));
1259 }
1260 std::sort(roots.begin(), roots.end(), vigra::detail::PolynomialRootCompare<Real>(eps));
1261 return true;
1262}
1263
1264template <class POLYNOMIAL, class VECTOR>
1265bool polynomialRootsEigenvalueMethod(POLYNOMIAL const & poly, VECTOR & roots)
1266{
1267 return polynomialRootsEigenvalueMethod(poly, roots, true);
1268}
1269
1270 /** Compute the real roots of a real polynomial using the eigenvalue method.
1271
1272 \a poly is a real polynomial (compatible to \ref vigra::PolynomialView),
1273 and \a roots a real valued vector (compatible to <tt>std::vector</tt>
1274 with a <tt>value_type</tt> compatible to the type <tt>POLYNOMIAL::Real</tt>) to which
1275 the roots are appended. The function calls \ref polynomialRootsEigenvalueMethod() and
1276 throws away all complex roots. It returns <tt>false</tt> if it fails to converge.
1277 The parameter <tt>polishRoots</tt> is ignored (it is only here for syntax compatibility
1278 with polynomialRealRoots()).
1279
1280
1281 <b>\#include</b> <vigra/eigensystem.hxx> or<br>
1282 <b>\#include</b> <vigra/linear_algebra.hxx><br>
1283 Namespaces: vigra and vigra::linalg
1284
1285 \see polynomialRealRoots(), vigra::Polynomial
1286 */
1287template <class POLYNOMIAL, class VECTOR>
1288bool polynomialRealRootsEigenvalueMethod(POLYNOMIAL const & p, VECTOR & roots, bool /* polishRoots */)
1289{
1290 typedef typename NumericTraits<typename VECTOR::value_type>::ComplexPromote Complex;
1291 ArrayVector<Complex> croots;
1292 if(!polynomialRootsEigenvalueMethod(p, croots))
1293 return false;
1294 for(unsigned int i = 0; i < croots.size(); ++i)
1295 if(croots[i].imag() == 0.0)
1296 roots.push_back(croots[i].real());
1297 return true;
1298}
1299
1300template <class POLYNOMIAL, class VECTOR>
1301bool polynomialRealRootsEigenvalueMethod(POLYNOMIAL const & p, VECTOR & roots)
1302{
1303 return polynomialRealRootsEigenvalueMethod(p, roots, true);
1304}
1305
1306
1307//@}
1308
1309} // namespace linalg
1310
1316
1317} // namespace vigra
1318
1319#endif // VIGRA_EIGENSYSTEM_HXX
size_type size() const
Definition array_vector.hxx:358
Definition array_vector.hxx:514
Base class for, and view to, MultiArray.
Definition multi_array.hxx:705
MultiArrayView< N-M, T, StrideTag > bindOuter(const TinyVector< Index, M > &d) const
Definition multi_array.hxx:2186
void copy(const MultiArrayView &rhs)
Definition multi_array.hxx:1218
Definition matrix.hxx:125
Linear algebra functions.
Definition eigensystem.hxx:50
Matrix< T, ALLLOC >::NormType norm(const Matrix< T, ALLLOC > &a)
bool symmetricEigensystem(MultiArrayView< 2, T, C1 > const &a, MultiArrayView< 2, T, C2 > &ew, MultiArrayView< 2, T, C3 > &ev)
Definition eigensystem.hxx:1008
MultiArrayIndex columnCount(const MultiArrayView< 2, T, C > &x)
Definition matrix.hxx:684
linalg::TemporaryMatrix< T > sq(MultiArrayView< 2, T, C > const &v)
bool nonsymmetricEigensystem(MultiArrayView< 2, T, C1 > const &a, MultiArrayView< 2, std::complex< T >, C2 > &ew, MultiArrayView< 2, T, C3 > &ev)
Definition eigensystem.hxx:1194
linalg::TemporaryMatrix< T > sqrt(MultiArrayView< 2, T, C > const &v)
linalg::TemporaryMatrix< T > abs(MultiArrayView< 2, T, C > const &v)
bool polynomialRootsEigenvalueMethod(POLYNOMIAL const &poly, VECTOR &roots, bool polishRoots)
Definition eigensystem.hxx:1233
MultiArrayIndex rowCount(const MultiArrayView< 2, T, C > &x)
Definition matrix.hxx:671
MultiArrayView< 2, T, C > columnVector(MultiArrayView< 2, T, C > const &m, MultiArrayIndex d)
Definition matrix.hxx:727
bool symmetricEigensystemNoniterative(MultiArrayView< 2, T, C1 > const &a, MultiArrayView< 2, T, C2 > &ew, MultiArrayView< 2, T, C3 > &ev)
Definition eigensystem.hxx:1152
bool polynomialRealRootsEigenvalueMethod(POLYNOMIAL const &p, VECTOR &roots, bool)
Definition eigensystem.hxx:1288
bool isSymmetric(const MultiArrayView< 2, T, C > &v)
Definition matrix.hxx:779
R imag(const FFTWComplex< R > &a)
imaginary part
Definition fftw3.hxx:1023
FFTWComplex< R >::NormType norm(const FFTWComplex< R > &a)
norm (= magnitude)
Definition fftw3.hxx:1037
R real(const FFTWComplex< R > &a)
real part
Definition fftw3.hxx:1016
void symmetric3x3Eigenvalues(T a00, T a01, T a02, T a11, T a12, T a22, T *r0, T *r1, T *r2)
Compute the eigenvalues of a 3x3 real symmetric matrix.
Definition mathutil.hxx:754
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.12.1 (Thu Feb 27 2025)