DGtal  1.5.beta
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  do
219  {
220  // Compute implicit shift
221  Quantity g = d[ l ];
222  Quantity p = ( d[ l + 1 ] - g ) / ( Quantity( 2.0 ) * e[ l ] );
223  Quantity r = Quantity( std::sqrt ( p * p + Quantity( 1.0 ) * Quantity( 1.0 )));
224  if( p < 0 ) r = -r;
225  d[ l ] = e[ l ] / ( p + r );
226  d[ l + 1 ] = e[ l ] * ( p + r );
227  Quantity dl1 = d[ l + 1 ];
228  Quantity h = g - d[ l ];
229  for( Dimension i = l + 2; i < dimension; ++i )
230  d[ i ] -= h;
231  f = f + h;
232 
233  // Implicit QL transformation.
234  p = d[ m ];
235  Quantity c = Quantity( 1.0 );
236  Quantity c2 = c;
237  Quantity c3 = c;
238  Quantity el1 = e[ l + 1 ];
239  Quantity s = Quantity( 0.0 );
240  Quantity s2 = Quantity( 0.0 );
241  for ( Dimension i = m - 1; i >= l && i <= m - 1; --i )
242  {
243  c3 = c2;
244  c2 = c;
245  s2 = s;
246  g = c * e[ i ];
247  h = c * p;
248  r = Quantity( std::sqrt ( p * p + e[ i ] * e[ i ] ));
249  e[ i + 1 ] = s * r;
250  s = e[ i ] / r;
251  c = p / r;
252  p = c * d[ i ] - s * g;
253  d[ i + 1 ] = h + s * ( c * g + s * d[ i ] );
254 
255  // Accumulate transformation.
256  for( Dimension k = 0; k < dimension; ++k )
257  {
258  h = V ( k, i + 1 );
259  V.setComponent ( k, i + 1, ( s * V ( k, i ) + c * h ));
260  V.setComponent ( k, i, ( c * V ( k, i ) - s * h ));
261  }
262  }
263 
264  p = - s * s2 * c3 * el1 * e[ l ] / dl1;
265  e[ l ] = s * p;
266  d[ l ] = c * p;
267  // Check for convergence.
268  }
269  while ( std::fabs ( e[ l ] ) > eps * tst1 );
270  }
271  d[ l ] = d[ l ] + f;
272  e[ l ] = Quantity( 0.0 );
273  }
274  // Sort eigenvalues and corresponding vectors.
275  for ( Dimension i = 0; i < dimensionMinusOne; ++i )
276  {
277  Dimension k = i;
278  Quantity p = d[ i ];
279 
280  for ( Dimension j = i + 1; j < dimension; ++j )
281  {
282  if ( d[ j ] < p )
283  {
284  k = j;
285  p = d[ j ];
286  }
287  }
288  if ( k != i )
289  {
290  d[ k ] = d[ i ];
291  d[ i ] = p;
292  for ( Dimension j = 0; j < dimension; ++j )
293  {
294  p = V ( j, i );
295  V.setComponent ( j, i, V ( j, k ));
296  V.setComponent ( j, k, p );
297  }
298  }
299  }
300 }
301 //-----------------------------------------------------------------------------
302 template <DGtal::Dimension TN, typename TComponent, typename TMatrix>
303 void
304 DGtal::EigenDecomposition<TN,TComponent,TMatrix>::
305 getEigenDecomposition( const Matrix& matrix, Matrix& eigenVectors, Vector& eigenValues )
306 {
307  Vector e; // Default constructor sets to zero vector;
308  eigenVectors = matrix; // copy matrix
309  eigenValues = e; // Sets to zero vector
310  tridiagonalize( eigenVectors, eigenValues, e );
311  decomposeQL( eigenVectors, eigenValues, e );
312 }
313 
314 // //
315 ///////////////////////////////////////////////////////////////////////////////
316 
317