DGtal  0.9.2
EigenDecomposition.ih
1 /**
2  * This program is free software: you can redistribute it and/or modify
3  * it under the terms of the GNU Lesser General Public License as
4  * published by the Free Software Foundation, either version 3 of the
5  * License, or (at your option) any later version.
6  *
7  * This program is distributed in the hope that it will be useful,
8  * but WITHOUT ANY WARRANTY; without even the implied warranty of
9  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
10  * GNU General Public License for more details.
11  *
12  * You should have received a copy of the GNU General Public License
13  * along with this program. If not, see <http://www.gnu.org/licenses/>.
14  *
15  **/
16 
17 /**
18  * @file EigenDecomposition.ih
19  * @author Jacques-Olivier Lachaud (\c jacques-olivier.lachaud@univ-savoie.fr )
20  * Laboratory of Mathematics (CNRS, UMR 5127), University of Savoie, France
21  *
22  * @date 2014/02/27
23  *
24  * Implementation of inline methods defined in EigenDecomposition.h
25  *
26  * This file is part of the DGtal library.
27  */
28 
29 
30 //////////////////////////////////////////////////////////////////////////////
31 #include <cstdlib>
32 //////////////////////////////////////////////////////////////////////////////
33 
34 ///////////////////////////////////////////////////////////////////////////////
35 // IMPLEMENTATION of inline methods.
36 ///////////////////////////////////////////////////////////////////////////////
37 
38 ///////////////////////////////////////////////////////////////////////////////
39 // ----------------------- Standard services ------------------------------
40 
41 //-----------------------------------------------------------------------------
42 template <DGtal::Dimension TN, typename TComponent, typename TMatrix>
43 void
44 DGtal::EigenDecomposition<TN,TComponent,TMatrix>::
45 tridiagonalize( Matrix& V, Vector& d, Vector& e )
46 {
47  for( Dimension j = 0; j < dimension; ++j )
48  d[ j ] = V ( dimensionMinusOne, j );
49  // Householder reduction to tridiagonal form.
50  for( Dimension i = dimensionMinusOne; i > 0 && i <= dimensionMinusOne; --i )
51  {
52  // Scale to avoid under/overflow.
53  Quantity scale = Quantity( 0.0 );
54  Quantity h = Quantity( 0.0 );
55  for( Dimension k = 0; k < i; ++k )
56  {
57  scale += Quantity( std::fabs( d[ k ] ));
58  }
59 
60  if( scale == Quantity( 0.0 ))
61  {
62  e[ i ] = d[ i - 1 ];
63 
64  for( Dimension j = 0; j < i; ++j )
65  {
66  d[ j ] = V(( i - 1 ), j );
67  V.setComponent ( i, j, Quantity( 0.0 ));
68  V.setComponent ( j, i, Quantity( 0.0 ));
69  }
70  }
71  else
72  {
73  // Generate Householder vector.
74  for ( Dimension k = 0; k < i; ++k )
75  {
76  d[ k ] /= scale;
77  h += d[ k ] * d[ k ];
78  }
79 
80  Quantity f = d[ i - 1 ];
81  Quantity g = Quantity( std::sqrt( h ));
82 
83  if ( f > Quantity( 0.0 ))
84  {
85  g = -g;
86  }
87 
88  e[ i ] = scale * g;
89  h -= f * g;
90  d[ i - 1 ] = f - g;
91 
92  for ( Dimension j = 0; j < i; ++j)
93  {
94  e[ j ] = Quantity( 0.0 );
95  }
96 
97  // Apply similarity transformation to remaining columns.
98  for ( Dimension j = 0; j < i; ++j )
99  {
100  f = d[ j ];
101  V.setComponent ( j, i, f );
102  g = e[ j ] + V( j, j ) * f;
103 
104  for ( Dimension k = j + 1; k <= i - 1; ++k )
105  {
106  g += V ( k, j ) * d[ k ];
107  e[ k ] += V ( k, j ) * f;
108  }
109 
110  e[ j ] = g;
111  }
112 
113  f = Quantity( 0.0 );
114  for ( Dimension j = 0; j < i; ++j )
115  {
116  e[ j ] /= h;
117  f += e[ j ] * d[ j ];
118  }
119 
120  double hh = f / ( h + h );
121 
122  for ( Dimension j = 0; j < i; ++j )
123  {
124  e[ j ] -= hh * d[ j ];
125  }
126 
127  for ( Dimension j = 0; j < i; ++j )
128  {
129  f = d[ j ];
130  g = e[ j ];
131 
132  for ( Dimension k = j; k <= i - 1; ++k )
133  {
134  V.setComponent ( k, j, V ( k, j ) - ( f * e[ k ] + g * d[ k ] ));
135  }
136 
137  d[ j ] = V( i - 1, j );
138  V.setComponent ( i, j, Quantity( 0.0 ));
139  }
140  }
141  d[ i ] = h;
142  }
143 
144  // Accumulate transformations.
145  for ( Dimension i = 0; i < dimensionMinusOne; ++i )
146  {
147  V.setComponent ( dimensionMinusOne, i, V ( i, i ));
148  V.setComponent ( i, i, Quantity( 1.0 ));
149  Quantity h = d[ i + 1 ];
150 
151  if ( h != Quantity( 0.0 ) )
152  {
153  for ( Dimension k = 0; k <= i; ++k )
154  {
155  d[ k ] = V ( k, i + 1 ) / h;
156  }
157 
158  for ( Dimension j = 0; j <= i; ++j )
159  {
160  Quantity g = Quantity( 0.0 );
161 
162  for ( Dimension k = 0; k <= i; ++k )
163  {
164  g += V ( k, i + 1 ) * V( k, j );
165  }
166 
167  for ( Dimension k = 0; k <= i; ++k )
168  {
169  V.setComponent ( k, j, V ( k, j ) - ( g * d[ k ] ));
170  }
171  }
172  }
173  for ( Dimension k = 0; k <= i; ++k )
174  {
175  V.setComponent ( k, i + 1, Quantity( 0.0 ));
176  }
177  }
178 
179  for ( Dimension j = 0; j < dimension; ++j )
180  {
181  d[ j ] = V ( dimensionMinusOne, j );
182  V.setComponent ( dimensionMinusOne, j, Quantity( 0.0 ));
183  }
184 
185  V.setComponent ( dimensionMinusOne, dimensionMinusOne, Quantity( 1.0 ));
186  e[ 0 ] = Quantity( 0.0 );
187 }
188 
189 //-----------------------------------------------------------------------------
190 template <DGtal::Dimension TN, typename TComponent, typename TMatrix>
191 void
192 DGtal::EigenDecomposition<TN,TComponent,TMatrix>::
193 decomposeQL( Matrix& V, Vector& d, Vector e )
194 {
195  for ( Dimension i = 1; i < dimension; ++i )
196  e[ i - 1 ] = e[ i ];
197 
198  e[ dimensionMinusOne ] = 0.0;
199 
200  Quantity f = Quantity( 0.0 );
201  Quantity tst1 = Quantity( 0.0 );
202  Quantity eps = Quantity( std::pow( 2.0, -52.0 ));
203  for( Dimension l = 0; l < dimension; ++l )
204  {
205  // Find small subdiagonal element
206  tst1 = Quantity( std::max( tst1, std::fabs ( d[ l ] ) + std::fabs( e[ l ] )));
207  Dimension m = l;
208  while ( m < dimension )
209  {
210  if ( std::fabs ( e[ m ] ) <= eps * tst1 ) break;
211  ++m;
212  }
213 
214  // If m == l, d[l] is an eigenvalue,
215  // otherwise, iterate.
216  if( m > l )
217  {
218  int iter = 0;
219  do
220  {
221  ++iter; // (Could check iteration count here.)
222  // Compute implicit shift
223  Quantity g = d[ l ];
224  Quantity p = ( d[ l + 1 ] - g ) / ( Quantity( 2.0 ) * e[ l ] );
225  Quantity r = Quantity( std::sqrt ( p * p + Quantity( 1.0 ) * Quantity( 1.0 )));
226  if( p < 0 ) r = -r;
227  d[ l ] = e[ l ] / ( p + r );
228  d[ l + 1 ] = e[ l ] * ( p + r );
229  Quantity dl1 = d[ l + 1 ];
230  Quantity h = g - d[ l ];
231  for( Dimension i = l + 2; i < dimension; ++i )
232  d[ i ] -= h;
233  f = f + h;
234 
235  // Implicit QL transformation.
236  p = d[ m ];
237  Quantity c = Quantity( 1.0 );
238  Quantity c2 = c;
239  Quantity c3 = c;
240  Quantity el1 = e[ l + 1 ];
241  Quantity s = Quantity( 0.0 );
242  Quantity s2 = Quantity( 0.0 );
243  for ( Dimension i = m - 1; i >= l && i <= m - 1; --i )
244  {
245  c3 = c2;
246  c2 = c;
247  s2 = s;
248  g = c * e[ i ];
249  h = c * p;
250  r = Quantity( std::sqrt ( p * p + e[ i ] * e[ i ] ));
251  e[ i + 1 ] = s * r;
252  s = e[ i ] / r;
253  c = p / r;
254  p = c * d[ i ] - s * g;
255  d[ i + 1 ] = h + s * ( c * g + s * d[ i ] );
256 
257  // Accumulate transformation.
258  for( Dimension k = 0; k < dimension; ++k )
259  {
260  h = V ( k, i + 1 );
261  V.setComponent ( k, i + 1, ( s * V ( k, i ) + c * h ));
262  V.setComponent ( k, i, ( c * V ( k, i ) - s * h ));
263  }
264  }
265 
266  p = - s * s2 * c3 * el1 * e[ l ] / dl1;
267  e[ l ] = s * p;
268  d[ l ] = c * p;
269  // Check for convergence.
270  }
271  while ( std::fabs ( e[ l ] ) > eps * tst1 );
272  }
273  d[ l ] = d[ l ] + f;
274  e[ l ] = Quantity( 0.0 );
275  }
276 
277  // Sort eigenvalues and corresponding vectors.
278  for ( Dimension i = 0; i < dimensionMinusOne; ++i )
279  {
280  Dimension k = i;
281  Quantity p = d[ i ];
282 
283  for ( Dimension j = i + 1; j < dimension; ++j )
284  {
285  if ( d[ j ] < p )
286  {
287  k = j;
288  p = d[ j ];
289  }
290  }
291  if ( k != i )
292  {
293  d[ k ] = d[ i ];
294  d[ i ] = p;
295  for ( Dimension j = 0; j < dimension; ++j )
296  {
297  p = V ( j, i );
298  V.setComponent ( j, i, V ( j, k ));
299  V.setComponent ( j, k, p );
300  }
301  }
302  }
303 }
304 //-----------------------------------------------------------------------------
305 template <DGtal::Dimension TN, typename TComponent, typename TMatrix>
306 void
307 DGtal::EigenDecomposition<TN,TComponent,TMatrix>::
308 getEigenDecomposition( const Matrix& matrix, Matrix& eigenVectors, Vector& eigenValues )
309 {
310  Vector e; // Default constructor sets to zero vector;
311  eigenVectors = matrix; // copy matrix
312  eigenValues = e; // Sets to zero vector
313  tridiagonalize( eigenVectors, eigenValues, e );
314  decomposeQL( eigenVectors, eigenValues, e );
315 }
316 
317 // //
318 ///////////////////////////////////////////////////////////////////////////////
319 
320