minimization.h
説明を見る。
1 //
2 // Copyright (c) 2003-2011, MIST Project, Nagoya University
3 // All rights reserved.
4 //
5 // Redistribution and use in source and binary forms, with or without modification,
6 // are permitted provided that the following conditions are met:
7 //
8 // 1. Redistributions of source code must retain the above copyright notice,
9 // this list of conditions and the following disclaimer.
10 //
11 // 2. Redistributions in binary form must reproduce the above copyright notice,
12 // this list of conditions and the following disclaimer in the documentation
13 // and/or other materials provided with the distribution.
14 //
15 // 3. Neither the name of the Nagoya University nor the names of its contributors
16 // may be used to endorse or promote products derived from this software
17 // without specific prior written permission.
18 //
19 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR
20 // IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
21 // FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
22 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
23 // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
24 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
25 // IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF
26 // THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
27 //
28 
38 
39 #ifndef __INCLUDE_MIST_MINIMIZATION__
40 #define __INCLUDE_MIST_MINIMIZATION__
41 
42 
43 #ifndef __INCLUDE_MIST_CONF_H__
44 #include "config/mist_conf.h"
45 #endif
46 
47 
48 #ifndef __INCLUDE_MIST_MATRIX__
49 #include "matrix.h"
50 #endif
51 
52 
53 #ifndef __INCLUDE_MIST_LIMITS__
54 #include "limits.h"
55 #endif
56 
57 
58 #include <cmath>
59 #include <vector>
60 #include <algorithm>
61 
62 // mist名前空間の始まり
63 _MIST_BEGIN
64 
65 
66 
67 namespace __minimization_utility__
68 {
69  template < class T, class Allocator, class Functor >
70  struct __convert_to_vector_functor__
71  {
72  typedef typename matrix< T, Allocator >::size_type size_type;
73  const matrix< T, Allocator > &ori_;
74  const matrix< T, Allocator > &dir_;
75  matrix< T, Allocator > &tmp_;
76  Functor &f_;
77 
78  __convert_to_vector_functor__( const matrix< T, Allocator > &ori, const matrix< T, Allocator > &dir, matrix< T, Allocator > &tmp, Functor &f ) : ori_( ori ), dir_( dir ), tmp_( tmp ), f_( f ){ }
79 
80  double operator ()( double x )
81  {
82  for( size_type i = 0 ; i < ori_.size( ) ; i++ )
83  {
84  tmp_[ i ] = ori_[ i ] + dir_[ i ] * x;
85  }
86  return( f_( tmp_ ) );
87  }
88  };
89 
90 
91  template < class Functor >
92  struct __no_copy_constructor_functor__
93  {
94  Functor &f_;
95  __no_copy_constructor_functor__( Functor &f ) : f_( f ){ }
96 
97  template < class PARAMETER >
98  double operator ()( const PARAMETER &x )
99  {
100  return( f_( x ) );
101  }
102 
103  template < class PARAMETER >
104  double operator ()( const PARAMETER &x ) const
105  {
106  return( f_( x ) );
107  }
108 
109  template < class PARAMETER >
110  double operator ()( const PARAMETER &x, size_t index )
111  {
112  return( f_( x, index ) );
113  }
114 
115  template < class PARAMETER >
116  double operator ()( const PARAMETER &x, size_t index ) const
117  {
118  return( f_( x, index ) );
119  }
120  };
121 
122 
123  template < class T, class Allocator, class Functor >
124  struct __gradient_vector_functor__
125  {
126  typedef typename matrix< T, Allocator >::size_type size_type;
127  typedef matrix< T, Allocator > matrix_type;
128  Functor &f_;
129  double d_;
130 
131  __gradient_vector_functor__( Functor &f, double d ) : f_( f ), d_( d ){ }
132 
133  const matrix_type operator ()( const matrix_type &v ) const
134  {
135  matrix_type dir( v.size( ), 1 ), tmp( v );
136  double len = 0.0, v1, v2;
137  size_type i;
138 
139  for( i = 0 ; i < dir.size( ) ; i++ )
140  {
141  tmp[ i ] = v[ i ] + d_;
142  v1 = f_( tmp );
143 
144  tmp[ i ] = v[ i ] - d_;
145  v2 = f_( tmp );
146 
147  tmp[ i ] = v[ i ];
148 
149  dir[ i ] = v1 - v2;
150  len += dir[ i ] * dir[ i ];
151  }
152 
153  if( len > 0 )
154  {
155  // 勾配方向ベクトルの正規化
156  len = std::sqrt( len );
157  for( i = 0 ; i < dir.size( ) ; i++ )
158  {
159  dir[ i ] /= len;
160  }
161  }
162 
163  return( dir );
164  }
165  };
166 
167  template < class T, class Allocator >
168  inline bool clipping( matrix< T, Allocator > &p1, matrix< T, Allocator > &p2, const double d1, const double d2, const double d )
169  {
170  typedef matrix< T, Allocator > matrix_type;
171 
172  bool c1 = d + d1 >= 0;
173  bool c2 = d + d2 >= 0;
174 
175  if( c1 && c2 )
176  {
177  return( true );
178  }
179  else if( !c1 && !c2 )
180  {
181  return( false );
182  }
183 
184  matrix_type v1 = p1;
185  matrix_type v2 = p2;
186  if( !c1 )
187  {
188  p1 = v2 + ( v1 - v2 ) * ( ( d + d2 ) / ( d2 - d1 ) );
189  }
190  if( !c2 )
191  {
192  p2 = v1 + ( v2 - v1 ) * ( ( d + d1 ) / ( d1 - d2 ) );
193  }
194 
195  return( true );
196  }
197 
198  template < class T, class Allocator >
199  inline bool clipping( matrix< T, Allocator > &p1, matrix< T, Allocator > &p2, const matrix< T, Allocator > &p, const matrix< T, Allocator > &dir, const matrix< T, Allocator > &box )
200  {
201  typedef matrix< T, Allocator > matrix_type;
202  typedef typename matrix_type::value_type value_type;
203  typedef typename matrix_type::size_type size_type;
204 
205  double infinity = 0.0;
206  size_type r;
207 
208  matrix_type offset( p.size( ), 1 );
209 
210  // もっとも長い対角方向の長さを取得し，バウンディングボックスの中央を求める
211  for( r = 0 ; r < box.rows( ) ; r++ )
212  {
213  double l = ( box( r, 0 ) - box( r, 1 ) );
214  infinity += l * l;
215  offset[ r ] = ( box( r, 0 ) + box( r, 1 ) ) / 2.0;
216  }
217 
218  infinity = std::sqrt( infinity );
219 
220  // まず，十分に遠い点を設定する
221  p1 = p - dir * infinity - offset;
222  p2 = p + dir * infinity - offset;
223 
224  bool flag = true;
225 
226  for( r = 0 ; r < box.rows( ) ; r++ )
227  {
228  // まず，下限をチェックする
229  flag = flag && clipping( p1, p2, p1[ r ], p2[ r ], std::abs( box( r, 0 ) - offset[ r ] ) );
230 
231  // 次に，上限をチェックする
232  flag = flag && clipping( p1, p2, -p1[ r ], -p2[ r ], std::abs( box( r, 1 ) - offset[ r ] ) );
233  }
234 
235  p1 += offset;
236  p2 += offset;
237 
238  return( flag );
239  }
240 
241  template < class T, class Allocator >
242  inline bool clipping_length( double &l1, double &l2, const matrix< T, Allocator > &p, const matrix< T, Allocator > &dir, const matrix< T, Allocator > &box )
243  {
244  typedef matrix< T, Allocator > matrix_type;
245  typedef typename matrix_type::value_type value_type;
246  typedef typename matrix_type::size_type size_type;
247 
248  matrix_type p1, p2;
249 
250  if( !clipping( p1, p2, p, dir, box ) )
251  {
252  return( false );
253  }
254 
255  l1 = l2 = 0.0;
256 
257  // もっとも長い対角方向の長さを取得し，バウンディングボックスの中央を求める
258  for( size_type r = 0 ; r < p1.size( ) ; r++ )
259  {
260  l1 += ( p[ r ] - p1[ r ] ) * ( p[ r ] - p1[ r ] );
261  l2 += ( p[ r ] - p2[ r ] ) * ( p[ r ] - p2[ r ] );
262  }
263 
264  // 数値計算誤差の関係で，バウンディングボックスを越えるのを回避する
265  l1 = - std::sqrt( l1 ) + 0.000000000001;
266  l2 = + std::sqrt( l2 ) - 0.000000000001;
267 
268  //std::cout << ( p + l1 * dir ).t( ) << std::endl;
269  //std::cout << ( p + l2 * dir ).t( ) << std::endl;
270 
271 
272  return( l1 < 0.0 || l2 > 0.0 );
273  }
274 
275 }
276 
277 
278 
286 
287 
288 
307 template < class Functor >
308 bool enclose( double &a, double &b, double &c, double &fa, double &fb, double &fc, Functor f, size_t max_iterations = 100 )
309 {
310  const double gold = ( 3.0 - std::sqrt( 5.0 ) ) / 2.0;
311  const double _1_gold = 1.0 / gold;
312  const double dust = type_limits< double >::tiny( ); // ゼロ除算を避けるためのゴミ値
313  const double limit = 100.0; // 放物線補外を適用する区間の最大更新幅
314 
315  if( a == b )
316  {
317  b = a + 1;
318  }
319 
320  fa = f( a );
321  fb = f( b );
322 
323  // fa < fb の場合は，aの方向に極小があると思われるので，aとbを入れ替える
324  if( fa < fb )
325  {
326  double tmp = a;
327  a = b;
328  b = tmp;
329  tmp = fa;
330  fa = fb;
331  fb = tmp;
332  }
333 
334  // 黄金分割比を用いて，次に探索する点cを決定する
335  c = a + _1_gold * ( b - a );
336  fc = f( c );
337 
338  // 区間の更新が不可能な場合の判定
339  if( fa == fb && fb == fc )
340  {
341  return( false );
342  }
343 
344  // f( a ) > f( b ) < f( c ) となるまで，区間の更新を続ける
345  // また，最大反復回数を超えた場合も終了する
346  size_t ite = 0;
347  while( fb > fc && ite++ < max_iterations )
348  {
349  // a, b, c, fa, fb, fc の値に，放物線を当てはめて極小が存在する位置を計算する
350  double ba = b - a;
351  double cb = c - b;
352  double fcb = fc - fb;
353  double fba = fb - fa;
354 
355  // ゼロ除算を防ぐために，ごみ値を足して，放物線補間を行う
356  double l1 = 2.0 * ( cb * fba - ba * fcb );
357  double l2 = std::abs( l1 );
358  double x = b + ( ba * ba * fcb + cb * cb * fba ) / ( l1 / l2 * ( l2 + dust ) );
359 
360  if( ( c - x ) * ( x - b ) > 0 )
361  {
362  // 補間点 x が区間 b < x < c を満たしていて，極小をはさむ区間を見つけることに成功
363  double fx = f( x );
364  if( fx < fc )
365  {
366  // 区間 ( b, c ) に極小が存在する
367  a = b;
368  b = x;
369  fa = fb;
370  fb = fx;
371  }
372  else
373  {
374  // 区間 ( a, x ) に極小が存在する
375  c = x;
376  fc = fx;
377  }
378 
379  // 区間を囲い込むことに成功したので終了する
380  break;
381  }
382  else if( ( b + limit * cb - x ) * ( x - c ) > 0 )
383  {
384  // 補間点 x が区間 b < c < x を満たしていて，許容範囲内に補間点が存在する
385  double fx = f( x );
386 
387  if( fx < fc )
388  {
389  // 現在の端である fc よりも関数値が小さい場合は，区間を更新する
390  // a <- b, b <- x てする
391  a = b;
392  b = x;
393  fa = fb;
394  fb = fx;
395 
396  // 新しい区間点 c を黄金分割比を用いて求める
397  c = a + _1_gold * ( b - a );
398  fc = f( c );
399  }
400  else
401  {
402  // b < c < x かつ fc < fx のため，区間の囲い込みに成功
403  // 区間 ( b, x ) に極小が存在する
404  a = b;
405  b = c;
406  c = x;
407  fa = fb;
408  fb = fc;
409  fc = fx;
410  }
411  }
412  else
413  {
414  // 補間点 x が区間 b < c < x を満たしているが，余りに補間点が遠すぎるので黄金分割する
415  // a <- b, b <- c として，新しい点 c を再計算する
416  a = b;
417  b = c;
418  c = a + _1_gold * ( b - a );
419  fa = fb;
420  fb = fc;
421  fc = f( c );
422  }
423  }
424 
425  return( ite < max_iterations );
426 }
427 
428 
430 namespace gold
431 {
451  template < class Functor >
452  double minimization( double a, double b, double &x, Functor f, double tolerance, size_t &iterations, size_t max_iterations, bool use_enclose = true )
453  {
454  double p, q, fp, fq;
455  const double gold = ( 3.0 - std::sqrt( 5.0 ) ) / 2.0;
456 
457  if( use_enclose )
458  {
459  p = b;
460  double fa, fb;
461  if( !enclose( a, p, b, fa, fp, fb, f ) )
462  {
463  return( fp );
464  }
465 
466  // a <= b <= c となるように区間を変更する
467  if( a > b )
468  {
469  double tmp = a;
470  a = b;
471  b = tmp;
472  tmp = fa;
473  fa = fb;
474  fb = tmp;
475  }
476 
477  // 区間の長いほうを決め，黄金分割によりもう一つの点を決定する
478  if( std::abs( p - a ) > std::abs( p - b ) )
479  {
480  q = ( p + a ) * 0.5;
481  fq = f( p );
482  }
483  else
484  {
485  q = ( p + b ) * 0.5;
486  fq = f( p );
487  }
488  }
489  else
490  {
491  if( a > b )
492  {
493  double tmp = a;
494  a = b;
495  b = tmp;
496  }
497  p = a + gold * ( b - a );
498  q = b - gold * ( b - a );
499 
500  fp = f( p );
501  fq = f( q );
502  }
503 
504  size_t ite = 0;
505  for( ; std::abs( a - b ) > tolerance * ( std::abs( p ) + std::abs( q ) ) && ite < max_iterations ; ite++ )
506  {
507  if( fp > fq )
508  {
509  // 区間 p < f( x ) < c の間に最小値が存在する
510  a = p;
511  p = q;
512  q = p + gold * ( b - p );
513  fp = fq;
514  fq = f( q );
515  }
516  else
517  {
518  // 区間 a < f( x ) < q の間に最小値が存在する
519  b = q;
520  q = p;
521  p = q - gold * ( q - a );
522  fq = fp;
523  fp = f( p );
524  }
525  }
526 
527  iterations = ite;
528 
529  if( fp < fq )
530  {
531  x = p;
532  return( fp );
533  }
534  else
535  {
536  x = q;
537  return( fq );
538  }
539  }
540 
541 
560  template < class Functor >
561  double minimization( double a, double b, double &x, Functor f, double tolerance, size_t max_iterations = 1000, bool use_enclose = true )
562  {
563  typedef __minimization_utility__::__no_copy_constructor_functor__< Functor > __no_copy_constructor_functor__;
564  size_t itenum = 0;
565  return( minimization( a, b, x, __no_copy_constructor_functor__( f ), tolerance, itenum, max_iterations, use_enclose ) );
566  }
567 }
568 
569 
570 
572 namespace brent
573 {
593  template < class Functor >
594  double minimization( double a, double b, double &x, Functor f, double tolerance, size_t &iterations, size_t max_iterations, bool use_enclose = true )
595  {
596  double u, v, w, fu, fv, fw, fx, e = 0.0, d = 0.0;
597  const double c = ( 3.0 - std::sqrt( 5.0 ) ) / 2.0;
598 
599  if( use_enclose )
600  {
601  double fa, fb;
602  x = b;
603  if( !enclose( a, x, b, fa, fx, fb, f ) )
604  {
605  return( fx );
606  }
607 
608  // a <= x <= b で極小を囲う区間に変更する
609  if( a > b )
610  {
611  double tmp = a;
612  a = b;
613  b = tmp;
614  tmp = fa;
615  fa = fb;
616  fb = tmp;
617  }
618 
619  v = a;
620  w = b;
621  u = x;
622  fv = fa;
623  fw = fb;
624  fu = fx;
625  }
626  else
627  {
628  // a <= x <= b で極小を囲う区間に変更する
629  if( a > b )
630  {
631  double tmp = a;
632  a = b;
633  b = tmp;
634  }
635 
636  v = w = x = a + c * ( b - a );
637  fv = fw = fx = f( x );
638  }
639 
640  size_t ite;
641  for( ite = 1 ; ite <= max_iterations ; ite++ )
642  {
643  double m = ( a + b ) * 0.5;
644  double tol = 1.0e-12 * std::abs( x ) + tolerance;
645  double t2 = 2.0 * tol;
646 
647  // 最後に判定した最小値候補点と，区間 [a, b] の中間との差が許容誤差内になったら終了する
648  if( std::abs( x - m ) <= t2 - 0.5 * ( b - a ) )
649  {
650  break;
651  }
652 
653  double p = 0.0, q = 0.0, r = 0.0;
654 
655  if( std::abs( e ) > tol )
656  {
657  // 放物線補間を行う
658  r = ( x - w ) * ( fx - fv );
659  q = ( x - v ) * ( fx - fw );
660  p = ( x - v ) * q - ( x - w ) * r;
661  q = 2.0 * ( q - r );
662 
663  if( q > 0 )
664  {
665  p = -p;
666  }
667  else
668  {
669  q = -q;
670  }
671 
672  r = e;
673  e = d;
674  }
675 
676  if( std::abs( p ) < std::abs( 0.5 * q * r ) && p > q * ( a - x ) && p < q * ( b - x ) )
677  //if( std::abs( p ) < std::abs( 0.5 * q * r ) && p < q * ( a - x ) && p < q * ( b - x ) )
678  {
679  // 放物線補間が適切に行われたので区間を更新する
680  d = p / q;
681  u = x + d;
682 
683  if( u - a < t2 || b - u < t2 )
684  {
685  d = x < m ? tol : -tol;
686  }
687  }
688  else
689  {
690  // 放物線補間は不適切なので黄金分割する
691  // 区間の大きいほうを黄金分割する
692  e = ( x < m ? b : a ) - x;
693  d = c * e;
694  }
695 
696  u = x + ( std::abs( d ) >= tol ? d : ( d > 0 ? tol : -tol ) );
697  fu = f( u );
698 
699  if( fu <= fx )
700  {
701  // より小さい値が見つかったので a, b, v, w を更新する
702  if( u < x )
703  {
704  // 区間 a < u < x に極小値がある
705  b = x;
706  }
707  else
708  {
709  // 区間 x < u < b に極小値がある
710  a = x;
711  }
712  v = w;
713  w = x;
714  fv = fw;
715  fw = fx;
716  x = u;
717  fx = fu;
718  }
719  else
720  {
721  // より大きい値が見つかったので a, b を更新する
722  if( u < x )
723  {
724  // 区間 u < x < b に極小値がある
725  a = u;
726  }
727  else
728  {
729  // 区間 a < x < u に極小値がある
730  b = u;
731  }
732 
733  // 新しく評価した点と従来評価した点の大小関係を調べる
734  if( fu <= fw || w == x )
735  {
736  v = w;
737  w = u;
738  fv = fw;
739  fw = fu;
740  }
741  else if( fu <= fv || v == x || v == w )
742  {
743  v = u;
744  fv = fu;
745  }
746  }
747  }
748 
749  iterations = ite;
750 
751  return( fx );
752  }
753 
754 
773  template < class Functor >
774  double minimization( double a, double b, double &x, Functor f, double tolerance, size_t max_iterations = 1000, bool use_enclose = true )
775  {
776  typedef __minimization_utility__::__no_copy_constructor_functor__< Functor > __no_copy_constructor_functor__;
777  size_t itenum = 0;
778  return( minimization( a, b, x, __no_copy_constructor_functor__( f ), tolerance, itenum, max_iterations, use_enclose ) );
779  }
780 }
781 
783 namespace Armijo
784 {
785  template < class T, class Allocator>
786  double mul( const array2< T, Allocator > &a, const array2< T, Allocator > &b )
787  {
788  double res = 0.0;
789  for ( size_t i = 0 ; i < a.size() ; i++ )
790  {
791  res += a[i] * b[i];
792  }
793  return res;
794  }
795 
796  template < class T, class Allocator>
797  double mul( const array2< rgb < T >, Allocator > &a, const array2< rgb < T >, Allocator > &b )
798  {
799  double res = 0.0;
800  for ( size_t i = 0 ; i < a.size() ; i++ )
801  {
802  res += a[i].r * b[i].r;
803  res += a[i].g * b[i].g;
804  res += a[i].b * b[i].b;
805  }
806  return res;
807  }
808 
822  template < class T, class Allocator,class Functor, class T2 >
823  double minimization( array2< T, Allocator > &p, Functor f, const array2< T, Allocator > &grad, const array2< T, Allocator > &d,T2 &data, const double rho = 0.7, const double c = 0.1, const size_t max_iterations = 10)
824  {
825  double alpha = 1.0;
826  double value = f( p, data );
827  double err = 0;
828  double fd = mul( grad, d );
829  array2< T, Allocator > start = p;
830 
831  int i = 0;
832 
833  do
834  {
835  if ( i == max_iterations )
836  {
837  p = start;
838  break;
839  }
840  if( i != 0 )
841  {
842  alpha *= rho;
843  }
844 
845  //p = start + alpha * d;
846  for( size_t k = 0 ; k < p.size() ; k++ )
847  {
848  p[k] = start[k] + alpha * d[k];
849  }
850 
851  i++;
852  //std::cout << "α = " << alpha ;
853  //std::cout << ", f = " << f( p, data ) << std::endl;
854  err = f( p, data );
855  }while( err > value + c * alpha * fd );
856 
857  return err;
858  }
859 
860 
873  template < class T, class Allocator,class Functor>
874  double minimization( array2< T, Allocator > &p, Functor f, const array2< T, Allocator > &grad, const array2< T, Allocator > &d, const double rho = 0.7, const double c = 0.1, const size_t max_iterations = 10)
875  {
876  double alpha = 1.0;
877  double value = f( p );
878  double err = 0;
879  double fd = mul( grad, d );
880  array2< T, Allocator > start = p;
881 
882  int i = 0;
883 
884  do
885  {
886  if ( i == max_iterations )
887  {
888  p = start;
889  break;
890  }
891  if( i != 0 )
892  {
893  alpha *= rho;
894  }
895 
896  //p = start + alpha * d;
897  for( size_t k = 0 ; k < p.size() ; k++ )
898  {
899  p[k] = start[k] + alpha * d[k];
900  }
901 
902  i++;
903  //std::cout << "α = " << alpha ;
904  //std::cout << ", f = " << f( p, data ) << std::endl;
905  err = f( p );
906  }while( err > value + c * alpha * fd );
907 
908  return err;
909  }
910 }
911 
913 namespace gradient_with_vector
914 {
928  template < class T, class Allocator, class Functor1, class Functor2 >
929  double minimization( matrix< T, Allocator > &p, Functor1 f, Functor2 g, double tolerance, size_t &iterations, size_t max_iterations = 1000 )
930  {
931  typedef typename matrix< T, Allocator >::value_type value_type;
932  typedef typename matrix< T, Allocator >::size_type size_type;
933  typedef typename matrix< T, Allocator >::difference_type difference_type;
934  typedef matrix< T, Allocator > matrix_type;
935 
936  matrix_type dir( p.size( ), 1 ), tmp( p.size( ), 1 );
937  double x, err, old_err = f( p );
938 
939  // 他変数関数を１変数関数に変換する
940  __minimization_utility__::__convert_to_vector_functor__< T, Allocator, Functor1 > functor( p, dir, tmp, f );
941 
942  size_t ite, i;
943  for( ite = 1 ; ite <= max_iterations ; ite++ )
944  {
945  // 勾配方向を計算する
946  double len = 0.0;
947  for( i = 0 ; i < dir.size( ) ; i++ )
948  {
949  // 勾配方向を計算する
950  dir[ i ] = g( p[ i ], i );
951  len += dir[ i ] * dir[ i ];
952  }
953 
954  if( len > 0 )
955  {
956  // 勾配方向ベクトルの正規化
957  len = std::sqrt( len );
958  for( i = 0 ; i < dir.size( ) ; i++ )
959  {
960  dir[ i ] /= len;
961  }
962  }
963  else
964  {
965  // 勾配の計算ができなくなったので終了する
966  break;
967  }
968 
969  // Brent の2次収束アルゴリズムを用いて dir 方向への最小化を行う
970  err = brent::minimization( 0.0, 1.0, x, functor, tolerance, 1000, true );
971 
972  //std::cout << p.t( ) << ", " << dir.t( ) << std::endl;
973 
974  if( 2.0 * std::abs( err - old_err ) <= tolerance * ( std::abs( err ) + std::abs( old_err ) ) || old_err <= err )
975  {
976  // 前回の最小化の結果からの変化量が、許容誤差内であったので終了する
977  if( err < old_err )
978  {
979  p += dir * x;
980  }
981  break;
982  }
983  else
984  {
985  old_err = err;
986  p += dir * x;
987  }
988  }
989 
990  iterations = ite;
991 
992  //std::cout << ite << std::endl;
993 
994  return( err );
995  }
996 
997 
1012  template < class T, class Allocator, class Functor1, class Functor2 >
1013  double minimization( matrix< T, Allocator > &p, const matrix< T, Allocator > &bound, Functor1 f, Functor2 g, double tolerance, size_t &iterations, size_t max_iterations = 1000 )
1014  {
1015  typedef typename matrix< T, Allocator >::value_type value_type;
1016  typedef typename matrix< T, Allocator >::size_type size_type;
1017  typedef typename matrix< T, Allocator >::difference_type difference_type;
1018  typedef matrix< T, Allocator > matrix_type;
1019 
1020  matrix_type dir( p.size( ), 1 ), tmp( p.size( ), 1 );
1021  double x, err, old_err = f( p );
1022 
1023  // 他変数関数を１変数関数に変換する
1024  __minimization_utility__::__convert_to_vector_functor__< T, Allocator, Functor1 > functor( p, dir, tmp, f );
1025 
1026  size_t ite, i;
1027  for( ite = 1 ; ite <= max_iterations ; ite++ )
1028  {
1029  // 勾配方向を計算する
1030  double len = 0.0;
1031  for( i = 0 ; i < dir.size( ) ; i++ )
1032  {
1033  // 勾配方向を計算する
1034  dir[ i ] = g( p[ i ], i );
1035  len += dir[ i ] * dir[ i ];
1036  }
1037 
1038  if( len > 0 )
1039  {
1040  // 勾配方向ベクトルの正規化
1041  len = std::sqrt( len );
1042  for( i = 0 ; i < dir.size( ) ; i++ )
1043  {
1044  dir[ i ] /= len;
1045  }
1046  }
1047  else
1048  {
1049  // 勾配の計算ができなくなったので終了する
1050  break;
1051  }
1052 
1053  double l1, l2;
1054  if( __minimization_utility__::clipping_length( l1, l2, p, dir, bound ) )
1055  {
1056  // Brent の2次収束アルゴリズムを用いて dir 方向への最小化を行う
1057  err = brent::minimization( l1, l2, x, functor, tolerance, 1000, false );
1058  }
1059 
1060  //std::cout << p.t( ) << ", " << dir.t( ) << std::endl;
1061 
1062 
1063  // 相対誤差を用いた収束判定
1064  if( 2.0 * std::abs( err - old_err ) <= tolerance * ( std::abs( err ) + std::abs( old_err ) ) || old_err <= err )
1065  {
1066  // 前回の最小化の結果からの変化量が、許容誤差内であったので終了する
1067  if( err < old_err )
1068  {
1069  p += dir * x;
1070  }
1071  break;
1072  }
1073  else
1074  {
1075  old_err = err;
1076  p += dir * x;
1077  }
1078  }
1079 
1080  iterations = ite;
1081 
1082  //std::cout << ite << std::endl;
1083 
1084  return( err );
1085  }
1086 
1087 
1101  template < class T, class Allocator, class Functor1, class Functor2 >
1102  double minimization( matrix< T, Allocator > &p, const matrix< T, Allocator > &bound, Functor1 f, Functor2 g, double tolerance, size_t max_iterations = 1000 )
1103  {
1104  typedef __minimization_utility__::__no_copy_constructor_functor__< Functor1 > __no_copy_constructor_functor1__;
1105  typedef __minimization_utility__::__no_copy_constructor_functor__< Functor2 > __no_copy_constructor_functor2__;
1106  size_t itenum = 0;
1107  return( minimization( p, bound, __no_copy_constructor_functor1__( f ), __no_copy_constructor_functor2__( g ), tolerance, itenum, max_iterations ) );
1108  }
1109 
1110 
1123  template < class T, class Allocator, class Functor1, class Functor2 >
1124  double minimization( matrix< T, Allocator > &p, Functor1 f, Functor2 g, double tolerance, size_t max_iterations = 1000 )
1125  {
1126  typedef __minimization_utility__::__no_copy_constructor_functor__< Functor1 > __no_copy_constructor_functor1__;
1127  typedef __minimization_utility__::__no_copy_constructor_functor__< Functor2 > __no_copy_constructor_functor2__;
1128  size_t itenum = 0;
1129  return( minimization( p, __no_copy_constructor_functor1__( f ), __no_copy_constructor_functor2__( g ), tolerance, itenum, max_iterations ) );
1130  }
1131 }
1132 
1133 
1135 namespace gradient
1136 {
1137 
1151  template < class T, class Allocator, class Functor >
1152  double minimization( matrix< T, Allocator > &p, Functor f, double tolerance, double distance, size_t &iterations, size_t max_iterations = 1000 )
1153  {
1154  typedef typename matrix< T, Allocator >::value_type value_type;
1155  typedef typename matrix< T, Allocator >::size_type size_type;
1156  typedef typename matrix< T, Allocator >::difference_type difference_type;
1157  typedef matrix< T, Allocator > matrix_type;
1158 
1159  matrix_type dir( p.size( ), 1 ), tmp = p;
1160  double x, v1, v2, err = 1.0e100, old_err = f( p );
1161  size_type i;
1162 
1163  // 他変数関数を１変数関数に変換する
1164  __minimization_utility__::__convert_to_vector_functor__< T, Allocator, Functor > functor( p, dir, tmp, f );
1165 
1166  size_t ite;
1167  for( ite = 1 ; ite <= max_iterations ; )
1168  {
1169  // 勾配方向を計算する
1170  double len = 0.0;
1171  for( i = 0 ; i < dir.size( ) ; i++ )
1172  {
1173  tmp[ i ] = p[ i ] + distance;
1174  v1 = f( tmp );
1175 
1176  tmp[ i ] = p[ i ] - distance;
1177  v2 = f( tmp );
1178 
1179  // 元に戻す
1180  tmp[ i ] = p[ i ];
1181 
1182  dir[ i ] = v2 - v1;
1183  len += dir[ i ] * dir[ i ];
1184  }
1185 
1186  if( len > 0 )
1187  {
1188  // 勾配方向ベクトルの正規化
1189  len = std::sqrt( len );
1190  for( i = 0 ; i < dir.size( ) ; i++ )
1191  {
1192  dir[ i ] /= len;
1193  }
1194  }
1195  else
1196  {
1197  // 勾配の計算ができなくなったので終了する
1198  break;
1199  }
1200 
1201  // Brent の2次収束アルゴリズムを用いて dir 方向への最小化を行う
1202  err = brent::minimization( 0.0, 1.0, x, functor, tolerance, 1000, true );
1203 
1204  //std::cout << err << ", " << p.t( ) << ", " << dir.t( ) << std::endl;
1205 
1206  ite++;
1207 
1208  if( ite > max_iterations || 2.0 * std::abs( old_err - err ) < tolerance * ( std::abs( old_err ) + std::abs( err ) ) || err >= old_err )
1209  {
1210  // 前回の最小化の結果からの変化量が、許容誤差内であったので終了する
1211  if( err < old_err )
1212  {
1213  p += dir * x;
1214  }
1215  break;
1216  }
1217  else
1218  {
1219  old_err = err;
1220  p += dir * x;
1221  }
1222  }
1223 
1224  iterations = ite;
1225  //std::cout << ite << std::endl;
1226 
1227  return( err );
1228  }
1229 
1230 
1245  template < class T, class Allocator, class Functor >
1246  double minimization( matrix< T, Allocator > &p, const matrix< T, Allocator > &bound, Functor f, double tolerance, double distance, size_t &iterations, size_t max_iterations = 1000 )
1247  {
1248  typedef typename matrix< T, Allocator >::value_type value_type;
1249  typedef typename matrix< T, Allocator >::size_type size_type;
1250  typedef typename matrix< T, Allocator >::difference_type difference_type;
1251  typedef matrix< T, Allocator > matrix_type;
1252 
1253  matrix_type dir( p.size( ), 1 ), tmp = p;
1254  double x = 0.0, v1, v2, err = 1.0e100, old_err = f( p );
1255  size_type i;
1256 
1257  // 他変数関数を１変数関数に変換する
1258  __minimization_utility__::__convert_to_vector_functor__< T, Allocator, Functor > functor( p, dir, tmp, f );
1259 
1260  size_t ite;
1261  for( ite = 1 ; ite <= max_iterations ; )
1262  {
1263  // 勾配方向を計算する
1264  double len = 0.0;
1265  for( i = 0 ; i < dir.size( ) ; i++ )
1266  {
1267  tmp[ i ] = p[ i ] + distance;
1268  v1 = f( tmp );
1269 
1270  tmp[ i ] = p[ i ] - distance;
1271  v2 = f( tmp );
1272 
1273  // 元に戻す
1274  tmp[ i ] = p[ i ];
1275 
1276  dir[ i ] = v2 - v1;
1277  len += dir[ i ] * dir[ i ];
1278  }
1279 
1280  if( len > 0 )
1281  {
1282  // 勾配方向ベクトルの正規化
1283  len = std::sqrt( len );
1284  for( i = 0 ; i < dir.size( ) ; i++ )
1285  {
1286  dir[ i ] /= len;
1287  }
1288  }
1289  else
1290  {
1291  // 勾配の計算ができなくなったので終了する
1292  break;
1293  }
1294 
1295  double l1 = 0.0, l2 = 0.0;
1296  if( __minimization_utility__::clipping_length( l1, l2, p, dir, bound ) )
1297  {
1298  // Brent の2次収束アルゴリズムを用いて dir 方向への最小化を行う
1299  err = brent::minimization( l1, l2, x, functor, tolerance, 1000, false );
1300  }
1301 
1302  //std::cout << err << ", " << p.t( ) << ", " << dir.t( ) << std::endl;
1303 
1304  ite++;
1305 
1306  if( ite > max_iterations || 2.0 * std::abs( old_err - err ) < tolerance * ( std::abs( old_err ) + std::abs( err ) ) )
1307  {
1308  // 前回の最小化の結果からの変化量が、許容誤差内であったので終了する
1309  if( err < old_err )
1310  {
1311  p += dir * x;
1312  }
1313  break;
1314  }
1315  else
1316  {
1317  old_err = err;
1318  p += dir * x;
1319  }
1320  }
1321 
1322  iterations = ite;
1323  //std::cout << ite << std::endl;
1324 
1325  return( err );
1326  }
1327 
1345  template < class T, class Allocator,class Functor1, class Functor2, class T2 >
1346  double minimization( array2< T, Allocator > &p, Functor1 f, Functor2 g, const T2 &data, double tolerance, size_t &iterations, const size_t max_iterations = 12, const double armijo_rho = 0.7, const double armijo_c = 0.1, const size_t armijo_max_iteration = 10 )
1347  {
1348  double beta = 0.0;
1349  double err = 0.0, old_err = f( p, data );
1350 
1351  array2< T, Allocator > grad( p.width(), p.height() ); // 勾配
1352 
1353  size_t ite;
1354  for( ite = 1 ; ite <= max_iterations ; ite++ )
1355  {
1356  //勾配方向を計算する
1357  g( p, grad, data );
1358 
1359  //Armijoの基準によりステップ幅α倍勾配方向へ画像を修正
1360  // x(k+1) = x(k) + α(k)d(k)
1361  err = Armijo::minimization( p, f, grad, -grad, data, armijo_rho, armijo_c, armijo_max_iteration );
1362 
1363  // 相対誤差を用いた収束判定
1364  if( 2.0 * std::abs( err - old_err ) <= tolerance * ( std::abs( err ) + std::abs( old_err ) ) || old_err <= err )
1365  {
1366  break;
1367  }
1368 
1369  old_err = err;
1370  std::cout << ite << " : Energy : " << err << std::endl;
1371  }
1372  iterations = ite;
1373  return err;
1374  }
1375 
1392  template < class T, class Allocator,class Functor1, class Functor2 >
1393  double minimization( array2< T, Allocator > &p, Functor1 f, Functor2 g, double tolerance, size_t &iterations, const size_t max_iterations = 12, const double armijo_rho = 0.7, const double armijo_c = 0.1, const size_t armijo_max_iteration = 10 )
1394  {
1395  double beta = 0.0;
1396  double err = 0.0, old_err = f( p );
1397 
1398  array2< T, Allocator > grad( p.width(), p.height() ); // 勾配
1399 
1400  size_t ite;
1401  for( ite = 1 ; ite <= max_iterations ; ite++ )
1402  {
1403  //勾配方向を計算する
1404  g( p, grad );
1405 
1406  //Armijoの基準によりステップ幅α倍勾配方向へ画像を修正
1407  // x(k+1) = x(k) + α(k)d(k)
1408  err = Armijo::minimization( p, f, grad, -grad, armijo_rho, armijo_c, armijo_max_iteration );
1409 
1410  // 相対誤差を用いた収束判定
1411  if( 2.0 * std::abs( err - old_err ) <= tolerance * ( std::abs( err ) + std::abs( old_err ) ) || old_err <= err )
1412  {
1413  break;
1414  }
1415 
1416  old_err = err;
1417  std::cout << ite << " : Energy : " << err << std::endl;
1418  }
1419  iterations = ite;
1420  return err;
1421  }
1434  template < class T, class Allocator, class Functor >
1435  double minimization( matrix< T, Allocator > &p, Functor f, double tolerance, double distance = 1.0, size_t max_iterations = 1000 )
1436  {
1437  typedef __minimization_utility__::__no_copy_constructor_functor__< Functor > __no_copy_constructor_functor__;
1438  size_t itenum = 0;
1439  return( minimization( p, __no_copy_constructor_functor__( f ), tolerance, distance, itenum, max_iterations ) );
1440  }
1441 
1442 
1456  template < class T, class Allocator, class Functor >
1457  double minimization( matrix< T, Allocator > &p, const matrix< T, Allocator > &bound, Functor f, double tolerance, double distance = 1.0, size_t max_iterations = 1000 )
1458  {
1459  typedef __minimization_utility__::__no_copy_constructor_functor__< Functor > __no_copy_constructor_functor__;
1460  size_t itenum = 0;
1461  return( minimization( p, bound, __no_copy_constructor_functor__( f ), tolerance, distance, itenum, max_iterations ) );
1462  }
1463 }
1464 
1465 
1467 namespace conjugate_gradient
1468 {
1469 
1470  template < class T, class Allocator >
1471  double norm( const array2< T, Allocator > a )
1472  {
1473  double res = 0.0;
1474  for( size_t i = 0 ; i < a.size() ; i++ )
1475  {
1476  res += a[i] * a[i];
1477  }
1478  return res;
1479  }
1480 
1481  template < class T, class Allocator >
1482  double norm( const array2< rgb< T >, Allocator > a )
1483  {
1484  double res = 0.0;
1485  for( size_t i = 0 ; i < a.size() ; i++ )
1486  {
1487  res += a[i].r * a[i].r;
1488  res += a[i].g * a[i].g;
1489  res += a[i].b * a[i].b;
1490  }
1491  return res;
1492  }
1493  template < class T, class Allocator>
1494  double mul( const array2< T, Allocator > &a, const array2< T, Allocator > &b )
1495  {
1496  double res = 0.0;
1497  for ( size_t i = 0 ; i < a.size() ; i++ )
1498  {
1499  res += a[i] * b[i];
1500  }
1501  return res;
1502  }
1503 
1504  template < class T, class Allocator>
1505  double mul( const array2< rgb < T >, Allocator > &a, const array2< rgb < T >, Allocator > &b )
1506  {
1507  double res = 0.0;
1508  for ( size_t i = 0 ; i < a.size() ; i++ )
1509  {
1510  res += a[i].r * b[i].r;
1511  res += a[i].g * b[i].g;
1512  res += a[i].b * b[i].b;
1513  }
1514  return res;
1515  }
1525  template < class T >
1526  double Polak_Ribiere_Polyak( const array2< T > &g0, const array2< T > &g1)
1527  {
1528  return mul( g1, g1 - g0 ) / pow( norm( g0 ), 2 );
1529  }
1538  template < class T >
1539  double Fletcher_Reeves( const array2< T > &g0, const array2< T > &g1)
1540  {
1541  return norm( g1 ) / norm( g0 );
1542  }
1543 
1544 
1561  template < class T, class Allocator,class Functor1, class Functor2>
1562  double minimization( array2< T, Allocator > &p, Functor1 f, Functor2 g, double tolerance, size_t &iterations, const size_t max_iterations = 12, const double armijo_rho = 0.7, const double armijo_c = 0.1, const size_t armijo_max_iteration = 10 )
1563  {
1564  double beta = 0.0;
1565  double err = 0.0, old_err = f( p );
1566 
1567  array2< T, Allocator > grad0( p.width(), p.height() ), grad1( p.width(), p.height() ); // 勾配
1568  array2< T, Allocator > d( p.width(), p.height() ); // 修正した勾配
1569 
1570  size_t ite;
1571  for( ite = 1 ; ite <= max_iterations ; ite++ )
1572  {
1573  //勾配方向を計算する
1574  g( p, grad1 );
1575 
1576  // 新しい勾配方向をβにより計算
1577  // p(k) = -g(k) + β(k)p(k-1)
1578  if( ite % p.width() == 1 )
1579  {
1580  d = - grad1;
1581  beta = 0.0;
1582  }
1583  else
1584  {
1585  // 修正した画像の勾配と修正する前の画像の勾配よりβを決定
1586  beta = Fletcher_Reeves( grad0, grad1 );
1587  //beta = Polak_Ribiere_Polyak( grad0, grad1 );
1588  //d = - grad1 + beta * d;
1589  for( size_t i = 0 ; i < d.size() ; i++ )
1590  {
1591  d[i] = -grad1[i] + beta * d[i];
1592  }
1593  }
1594 
1595  //Armijoの基準によりステップ幅α倍勾配方向へ画像を修正
1596  err = Armijo::minimization( p, f, grad1, d, armijo_rho, armijo_c, armijo_max_iteration );
1597 
1598  // 相対誤差を用いた収束判定
1599  if( 2.0 * std::abs( err - old_err ) <= tolerance * ( std::abs( err ) + std::abs( old_err ) ) || old_err <= err )
1600  {
1601  break;
1602  }
1603 
1604  old_err = err;
1605  grad0 = grad1;
1606  std::cout << ite << " : Energy : " << err << " : β -> " << beta << std::endl;
1607  }
1608  iterations = ite;
1609  return err;
1610  }
1611 
1629  template < class T, class Allocator,class Functor1, class Functor2, class T2 >
1630  double minimization( array2< T, Allocator > &p, Functor1 f, Functor2 g, const T2 &data, double tolerance, size_t &iterations, const size_t max_iterations = 12, const double armijo_rho = 0.7, const double armijo_c = 0.1, const size_t armijo_max_iteration = 10 )
1631  {
1632  double beta = 0.0;
1633  double err = 0.0, old_err = f( p, data );
1634 
1635  array2< T, Allocator > grad0( p.width(), p.height() ), grad1( p.width(), p.height() ); // 勾配
1636  array2< T, Allocator > d( p.width(), p.height() ); // 修正した勾配
1637 
1638  size_t ite;
1639  for( ite = 1 ; ite <= max_iterations ; ite++ )
1640  {
1641  //勾配方向を計算する
1642  g( p, grad1, data );
1643 
1644  // 新しい勾配方向をβにより計算
1645  // p(k) = -g(k) + β(k)p(k-1)
1646  if( ite % p.width() == 1 )
1647  {
1648  d = - grad1;
1649  beta = 0.0;
1650  }
1651  else
1652  {
1653  // 修正した画像の勾配と修正する前の画像の勾配よりβを決定
1654  beta = Fletcher_Reeves( grad0, grad1 );
1655  //beta = Polak_Ribiere_Polyak( grad0, grad1 );
1656  //d = - grad1 + beta * d;
1657  for( size_t i = 0 ; i < d.size() ; i++ )
1658  {
1659  d[i] = -grad1[i] + beta * d[i];
1660  }
1661  }
1662 
1663  //Armijoの基準によりステップ幅α倍勾配方向へ画像を修正
1664  err = Armijo::minimization( p, f, grad1, d, data, armijo_rho, armijo_c, armijo_max_iteration );
1665 
1666  // 相対誤差を用いた収束判定
1667  if( 2.0 * std::abs( err - old_err ) <= tolerance * ( std::abs( err ) + std::abs( old_err ) ) || old_err <= err )
1668  {
1669  break;
1670  }
1671 
1672  old_err = err;
1673  grad0 = grad1;
1674  std::cout << ite << " : Energy : " << err << " : β -> " << beta << std::endl;
1675  }
1676  iterations = ite;
1677  return err;
1678  }
1679 
1680 }
1682 namespace powell
1683 {
1697  template < class T, class Allocator, class Functor >
1698  double minimization( matrix< T, Allocator > &p, matrix< T, Allocator > &dirs, Functor f, double tolerance, size_t &iterations, size_t max_iterations = 1000 )
1699  {
1700  typedef typename matrix< T, Allocator >::value_type value_type;
1701  typedef typename matrix< T, Allocator >::size_type size_type;
1702  typedef typename matrix< T, Allocator >::difference_type difference_type;
1703  typedef matrix< T, Allocator > matrix_type;
1704 
1705  matrix_type odirs( dirs ), dir( p.size( ), 1 ), tmp( p.size( ), 1 ), p0( p ), pn( p );
1706  double x, fp = f( p ), fp0 = fp, delta;
1707 
1708  // 他変数関数を１変数関数に変換する
1709  __minimization_utility__::__convert_to_vector_functor__< T, Allocator, Functor > functor( p, dir, tmp, f );
1710 
1711  size_t ite;
1712  size_type r, c;
1713  for( ite = 0 ; ite < max_iterations ; ite++ )
1714  {
1715  // 探索開始前の評価値を覚えておく
1716  fp0 = fp;
1717  delta = -1.0;
1718  size_type index = 0;
1719 
1720  for( c = 0 ; c < dirs.cols( ) ; c++ )
1721  {
1722  // 探索に用いる方向集合をコピーする
1723  for( r = 0 ; r < dirs.rows( ) ; r++ )
1724  {
1725  dir[ r ] = dirs( r, c );
1726  }
1727 
1728  double old_fp = fp;
1729 
1730  // Brent の2次収束アルゴリズムを用いて dir 方向への最小化を行う
1731  fp = brent::minimization( 0.0, 1.0, x, functor, tolerance, 1000, true );
1732 
1733  for( r = 0 ; r < p.size( ) ; r++ )
1734  {
1735  p[ r ] += dir[ r ] * x;
1736  }
1737 
1738  double d = std::abs( fp - old_fp );
1739  if( d > delta )
1740  {
1741  index = c;
1742  delta = d;
1743  }
1744  }
1745 
1746  // 相対誤差を用いた収束判定
1747  if( 2.0 * std::abs( fp - fp0 ) <= tolerance * ( std::abs( fp ) + std::abs( fp0 ) ) || fp0 <= fp )
1748  {
1749  break;
1750  }
1751 
1752  // Acton の方法を用いて，新しい方向集合を求める
1753  if( ite < max_iterations )
1754  {
1755  // 新しい方向を求める
1756  for( r = 0 ; r < p.size( ) ; r++ )
1757  {
1758  dir[ r ] = p[ r ] - p0[ r ];
1759  pn[ r ] = 2.0 * p[ r ] - p0[ r ];
1760  }
1761 
1762  double fe = f( pn );
1763 
1764  if( fe < fp0 )
1765  {
1766  // 現在の方向集合を更新する
1767  double tmp = fp0 - fp - delta;
1768  double ttt = 2.0 * ( fp0 - 2.0 * fp + fe ) * tmp * tmp - delta * ( fp0 - fe ) * ( fp0 - fe );
1769  if( ttt < 0 )
1770  {
1771  // Brent の2次収束アルゴリズムを用いて，新しい dir 方向への最小化を行う
1772  fp = brent::minimization( 0.0, 1.0, x, functor, tolerance, 1000, true );
1773 
1774  // 現在のパラメータ更新
1775  for( r = 0 ; r < p.size( ) ; r++ )
1776  {
1777  p[ r ] += dir[ r ] * x;
1778  }
1779 
1780  // 方向集合の一番最後に，新しい方向を追加する
1781  if( index < dirs.rows( ) - 1 )
1782  {
1783  for( r = 0 ; r < dirs.rows( ) ; r++ )
1784  {
1785  dirs( r, index ) = dirs( r, dirs.rows( ) - 1 );
1786  dirs( r, dirs.rows( ) - 1 ) = dir[ r ];
1787  }
1788  }
1789  else
1790  {
1791  for( r = 0 ; r < dirs.rows( ) ; r++ )
1792  {
1793  dirs( r, index ) = dir[ r ];
1794  }
1795  }
1796  }
1797  }
1798 
1799  // 新しい方向を求める
1800  for( r = 0 ; r < p.size( ) ; r++ )
1801  {
1802  p0[ r ] = p[ r ];
1803  }
1804  }
1805  }
1806 
1807  iterations = ite;
1808 
1809  return( fp );
1810  }
1811 
1812 
1827  template < class T, class Allocator, class Functor >
1828  double minimization( matrix< T, Allocator > &p, matrix< T, Allocator > &dirs, matrix< T, Allocator > &bound, Functor f, double tolerance, size_t &iterations, size_t max_iterations = 1000 )
1829  {
1830  typedef typename matrix< T, Allocator >::value_type value_type;
1831  typedef typename matrix< T, Allocator >::size_type size_type;
1832  typedef typename matrix< T, Allocator >::difference_type difference_type;
1833  typedef matrix< T, Allocator > matrix_type;
1834 
1835  matrix_type dir( p.size( ), 1 ), tmp( p.size( ), 1 ), p0( p ), pn( p );
1836  double x = 0.0, fp0 = 1.0e100, fp = f( p ), delta;
1837  double l1 = 0.0, l2 = 0.0;
1838 
1839  // 他変数関数を１変数関数に変換する
1840  __minimization_utility__::__convert_to_vector_functor__< T, Allocator, Functor > functor( p, dir, tmp, f );
1841 
1842  size_t ite;
1843  size_type r, c;
1844  for( ite = 1 ; ite <= max_iterations ; ite++ )
1845  {
1846  // 探索開始前の評価値を覚えておく
1847  fp0 = fp;
1848  delta = 0.0;
1849  size_type index = 0;
1850 
1851  for( c = 0 ; c < dirs.cols( ) ; c++ )
1852  {
1853  // 探索に用いる方向集合をコピーする
1854  for( r = 0 ; r < dirs.rows( ) ; r++ )
1855  {
1856  dir[ r ] = dirs( r, c );
1857  }
1858 
1859  double old_fp = fp;
1860 
1861 
1862  if( __minimization_utility__::clipping_length( l1, l2, p, dir, bound ) )
1863  {
1864  // Brent の2次収束アルゴリズムを用いて dir 方向への最小化を行う
1865  double nfp = brent::minimization( l1, l2, x, functor, tolerance, 1000, false );
1866 
1867  // 関数値が減少した場合のみ値を更新する
1868  if( nfp < fp )
1869  {
1870  fp = nfp;
1871 
1872  for( r = 0 ; r < p.size( ) ; r++ )
1873  {
1874  p[ r ] += dir[ r ] * x;
1875  }
1876 
1877  double d = std::abs( fp - old_fp );
1878  if( d > delta )
1879  {
1880  index = c;
1881  delta = d;
1882  }
1883  }
1884  }
1885  }
1886 
1887  // 相対誤差を用いた収束判定
1888  if( 2.0 * std::abs( fp - fp0 ) <= tolerance * ( std::abs( fp ) + std::abs( fp0 ) ) )
1889  {
1890  break;
1891  }
1892 
1893  // Acton の方法を用いて，新しい方向集合を求める
1894  if( ite <= max_iterations )
1895  {
1896  // 新しい方向を求める
1897  for( r = 0 ; r < p.size( ) ; r++ )
1898  {
1899  dir[ r ] = p[ r ] - p0[ r ];
1900  pn[ r ] = 2.0 * p[ r ] - p0[ r ];
1901  }
1902 
1903  double fe = f( pn );
1904 
1905  if( fe < fp0 )
1906  {
1907  // 現在の方向集合を更新する
1908  double tmp = fp0 - fp - delta;
1909  double ttt = 2.0 * ( fp0 - 2.0 * fp + fe ) * tmp * tmp - delta * ( fp0 - fe ) * ( fp0 - fe );
1910  if( ttt < 0 )
1911  {
1912  if( __minimization_utility__::clipping_length( l1, l2, p, dir, bound ) )
1913  {
1914  // Brent の2次収束アルゴリズムを用いて，新しい dir 方向への最小化を行う
1915  double nfp = brent::minimization( l1, l2, x, functor, tolerance, 1000, false );
1916 
1917  // 関数値が減少した場合のみ値を更新する
1918  if( nfp < fp )
1919  {
1920  fp = nfp;
1921  for( r = 0 ; r < p.size( ) ; r++ )
1922  {
1923  p[ r ] += dir[ r ] * x;
1924  }
1925  }
1926  }
1927 
1928  // 方向集合の一番最後に，新しい方向を追加する
1929  if( index < dirs.rows( ) - 1 )
1930  {
1931  for( r = 0 ; r < dirs.rows( ) ; r++ )
1932  {
1933  dirs( r, index ) = dirs( r, dirs.rows( ) - 1 );
1934  dirs( r, dirs.rows( ) - 1 ) = dir[ r ];
1935  }
1936  }
1937  else
1938  {
1939  for( r = 0 ; r < dirs.rows( ) ; r++ )
1940  {
1941  dirs( r, index ) = dir[ r ];
1942  }
1943  }
1944  }
1945  }
1946 
1947  // 新しい方向を求める
1948  for( r = 0 ; r < p.size( ) ; r++ )
1949  {
1950  p0[ r ] = p[ r ];
1951  }
1952  }
1953  }
1954 
1955  iterations = ite;
1956 
1957  return( fp );
1958  }
1959 
1972  template < class T, class Allocator, class Functor >
1973  double minimization( matrix< T, Allocator > &p, matrix< T, Allocator > &dirs, Functor f, double tolerance, size_t max_iterations = 1000 )
1974  {
1975  typedef __minimization_utility__::__no_copy_constructor_functor__< Functor > __no_copy_constructor_functor__;
1976  size_t itenum = 0;
1977  return( minimization( p, dirs, __no_copy_constructor_functor__( f ), tolerance, itenum, max_iterations ) );
1978  }
1979 
1993  template < class T, class Allocator, class Functor >
1994  double minimization( matrix< T, Allocator > &p, matrix< T, Allocator > &dirs, matrix< T, Allocator > &bound, Functor f, double tolerance, size_t max_iterations = 1000 )
1995  {
1996  typedef __minimization_utility__::__no_copy_constructor_functor__< Functor > __no_copy_constructor_functor__;
1997  size_t itenum = 0;
1998  return( minimization( p, dirs, bound, __no_copy_constructor_functor__( f ), tolerance, itenum, max_iterations ) );
1999  }
2000 }
2001 
2002 
2004 namespace lucidi
2005 {
2006  template < class Functor >
2007  inline double expantion_step( double alpha, double alpha_max, double delta, double gamma, double fp0, Functor f )
2008  {
2009  while( true )
2010  {
2011  if( alpha >= alpha_max )
2012  {
2013  return( alpha_max );
2014  }
2015  else
2016  {
2017  double a = alpha / delta;
2018  if( f( a ) > fp0 - gamma * a * a )
2019  {
2020  return( alpha );
2021  }
2022 
2023  alpha = a;
2024  }
2025  }
2026  }
2027 
2028 
2042  template < class T, class Allocator, class Functor >
2043  double minimization( matrix< T, Allocator > &p, matrix< T, Allocator > &dirs, Functor f, double tolerance, size_t &iterations, size_t max_iterations = 1000 )
2044  {
2045  typedef typename matrix< T, Allocator >::value_type value_type;
2046  typedef typename matrix< T, Allocator >::size_type size_type;
2047  typedef typename matrix< T, Allocator >::difference_type difference_type;
2048  typedef matrix< T, Allocator > matrix_type;
2049 
2050  matrix_type dir( p.size( ), 1 ), tmp( p.size( ), 1 ), p0( p ), a( p.size( ), 1 );
2051  double fp0 = 1.0e100, fp = f( p );
2052 
2053  // 他変数関数を１変数関数に変換する
2054  __minimization_utility__::__convert_to_vector_functor__< T, Allocator, Functor > functor( p, dir, tmp, f );
2055 
2056  const double __alpha__ = 0.5;
2057  const double theta = 0.5;
2058  const double gamma = 1.0e-6;
2059  const double delta = 0.25;
2060  const double infinity = type_limits< double >::maximum( );
2061 
2062  size_t ite;
2063  size_type r, c;
2064 
2065  for( r = 0 ; r < a.size( ) ; r++ )
2066  {
2067  a[ r ] = __alpha__;
2068  }
2069 
2070  for( ite = 1 ; ite <= max_iterations ; ite++ )
2071  {
2072  // 探索開始前の評価値を覚えておく
2073  fp0 = fp;
2074 
2075  size_type alpha_count = 0;
2076 
2077  for( c = 0 ; c < dirs.cols( ) ; c++ )
2078  {
2079  double alpha = a[ c ];
2080 
2081  // 探索に用いる方向集合をコピーする
2082  for( r = 0 ; r < dirs.rows( ) ; r++ )
2083  {
2084  dir[ r ] = dirs( r, c );
2085  }
2086 
2087  if( alpha > 0.0 )
2088  {
2089  double gaa = gamma * alpha * alpha;
2090  if( functor( alpha ) <= fp0 - gaa )
2091  {
2092  a[ c ] = alpha = expantion_step( alpha, infinity, delta, gamma, fp0, functor );
2093 
2094  // 位置を更新
2095  for( r = 0 ; r < p.size( ) ; r++ )
2096  {
2097  p[ r ] += dir[ r ] * alpha;
2098  }
2099  }
2100  else if( functor( -alpha ) <= fp0 - gaa )
2101  {
2102  // 探索に用いる方向集合をコピーする
2103  for( r = 0 ; r < dirs.rows( ) ; r++ )
2104  {
2105  dir[ r ] = dirs( r, c ) = -dir[ r ];
2106  }
2107 
2108  a[ c ] = alpha = expantion_step( alpha, infinity, delta, gamma, fp0, functor );
2109 
2110  // 位置を更新
2111  for( r = 0 ; r < p.size( ) ; r++ )
2112  {
2113  p[ r ] += dir[ r ] * alpha;
2114  }
2115  }
2116  else
2117  {
2118  a[ c ] = theta * alpha;
2119  alpha = 0.0;
2120  alpha_count++;
2121  }
2122  }
2123  else
2124  {
2125  a[ c ] = theta * alpha;
2126  alpha = 0.0;
2127  alpha_count++;
2128  }
2129  }
2130 
2131  fp = f( p );
2132 
2133  if( alpha_count < dirs.cols( ) )
2134  {
2135  // 相対誤差を用いた収束判定
2136  if( 2.0 * std::abs( fp - fp0 ) <= tolerance * ( std::abs( fp ) + std::abs( fp0 ) ) )
2137  {
2138  break;
2139  }
2140  }
2141  }
2142 
2143  iterations = ite;
2144 
2145  return( fp );
2146  }
2147 
2148 
2161  template < class T, class Allocator, class Functor >
2162  double minimization( matrix< T, Allocator > &p, matrix< T, Allocator > &dirs, Functor f, double tolerance, size_t max_iterations = 1000 )
2163  {
2164  typedef __minimization_utility__::__no_copy_constructor_functor__< Functor > __no_copy_constructor_functor__;
2165  size_t itenum = 0;
2166  return( minimization( p, dirs, __no_copy_constructor_functor__( f ), tolerance, itenum, max_iterations ) );
2167  }
2168 }
2169 
2170 
2171 
2172 
2178 namespace condor
2179 {
2180  // CONDORアルゴリズム内部で利用する関数など
2181  namespace __condor_utility__
2182  {
2183  inline double minimum( double v1, double v2 )
2184  {
2185  return( v1 < v2 ? v1 : v2 );
2186  }
2187 
2188  inline double maximum( double v1, double v2 )
2189  {
2190  return( v1 > v2 ? v1 : v2 );
2191  }
2192 
2193  inline double minimum( double v1, double v2, double v3 )
2194  {
2195  if( v1 < v2 )
2196  {
2197  if( v1 < v3 )
2198  {
2199  return( v1 );
2200  }
2201  else
2202  {
2203  return( v3 );
2204  }
2205  }
2206  else if( v2 < v3 )
2207  {
2208  return( v2 );
2209  }
2210  else
2211  {
2212  return( v3 );
2213  }
2214  }
2215 
2216  inline double minimum( double v1, double v2, double v3, double v4 )
2217  {
2218  return( minimum( minimum( v1, v2 ), minimum( v3, v4 ) ) );
2219  }
2220 
2221  inline double maximum( double v1, double v2, double v3 )
2222  {
2223  if( v1 > v2 )
2224  {
2225  if( v1 > v3 )
2226  {
2227  return( v1 );
2228  }
2229  else
2230  {
2231  return( v3 );
2232  }
2233  }
2234  else if( v2 > v3 )
2235  {
2236  return( v2 );
2237  }
2238  else
2239  {
2240  return( v3 );
2241  }
2242  }
2243 
2244  inline double maximum( double v1, double v2, double v3, double v4 )
2245  {
2246  return( maximum( maximum( v1, v2 ), maximum( v3, v4 ) ) );
2247  }
2248 
2249  template < class Matrix >
2250  inline double frobenius_norm( const Matrix &H )
2251  {
2252  typedef Matrix matrix_type;
2253  typedef typename matrix_type::size_type size_type;
2254 
2255  double val = 0.0;
2256  for( size_type r = 0 ; r < H.size( ) ; r++ )
2257  {
2258  val += H[ r ] * H[ r ];
2259  }
2260 
2261  return( std::sqrt( val ) );
2262  }
2263 
2264  template < class Matrix >
2265  inline double infinitum_norm( const Matrix &H )
2266  {
2267  typedef Matrix matrix_type;
2268  typedef typename matrix_type::size_type size_type;
2269 
2270  double max = 0.0;
2271  for( size_type c = 0 ; c < H.cols( ) ; c++ )
2272  {
2273  double val = 0.0;
2274  for( size_type r = 0 ; r < H.rows( ) ; r++ )
2275  {
2276  val += std::abs( H( r, c ) );
2277  }
2278 
2279  if( max < val )
2280  {
2281  max = val;
2282  }
2283  }
2284 
2285  return( max );
2286  }
2287 
2288  template < class Matrix >
2289  inline double inner_product( const Matrix &m1, const Matrix &m2 )
2290  {
2291  typedef Matrix matrix_type;
2292  typedef typename matrix_type::value_type value_type;
2293  typedef typename matrix_type::size_type size_type;
2294  typedef typename matrix_type::difference_type difference_type;
2295 
2296  double sum = 0.0;
2297  for( size_type i = 0 ; i < m1.size( ) ; i++ )
2298  {
2299  sum += m1[ i ] * m2[ i ];
2300  }
2301 
2302  return( sum );
2303  }
2304 
2305  template < class Matrix >
2306  inline double inner_product( const Matrix &m1, const Matrix &H, const Matrix &m2 )
2307  {
2308  typedef Matrix matrix_type;
2309  typedef typename matrix_type::size_type size_type;
2310 
2311  double sum = 0.0;
2312  for( size_type c = 0 ; c < H.cols( ) ; c++ )
2313  {
2314  double val = 0.0;
2315  for( size_type r = 0 ; r < H.rows( ) ; r++ )
2316  {
2317  val += m1[ r ] * H( r, c );
2318  }
2319 
2320  sum += val * m2[ c ];
2321  }
2322  return( sum );
2323  }
2324 
2325  template < class Matrix >
2326  inline double inner_product( const Matrix &m1, const Matrix &H, const Matrix &m2, double lambda )
2327  {
2328  typedef Matrix matrix_type;
2329  typedef typename matrix_type::size_type size_type;
2330 
2331  double sum = 0.0;
2332  for( size_type c = 0 ; c < H.cols( ) ; c++ )
2333  {
2334  double val = lambda * m1[ c ];
2335  for( size_type r = 0 ; r < H.rows( ) ; r++ )
2336  {
2337  val += m1[ r ] * H( r, c );
2338  }
2339 
2340  sum += val * m2[ c ];
2341  }
2342  return( sum );
2343  }
2344 
2345  struct __index_value_pair__
2346  {
2347  size_t index;
2348  double value;
2349 
2350  __index_value_pair__( size_t indx, double val ) : index( indx ), value( val ){ }
2351 
2352  bool operator <( const __index_value_pair__ &v ) const
2353  {
2354  // 距離値の大きい順に並ぶようにする
2355  return( value > v.value );
2356  }
2357  };
2358 
2359  template < class Matrix >
2360  inline void solve( const Matrix &A, Matrix &b )
2361  {
2362  typedef Matrix matrix_type;
2363  typedef typename matrix_type::size_type size_type;
2364  typedef typename matrix_type::difference_type difference_type;
2365 
2366  for( size_type r = 0 ; r < A.rows( ) ; r++ )
2367  {
2368  double sum = b[ r ];
2369  for( difference_type c = r - 1 ; c >= 0 ; c-- )
2370  {
2371  sum -= A( r, c ) * b[ c ];
2372  }
2373 
2374  b[ r ] = sum / A( r, r );
2375  }
2376  }
2377 
2378  template < class Matrix >
2379  inline void solve_( const Matrix &A, Matrix &b )
2380  {
2381  typedef Matrix matrix_type;
2382  typedef typename matrix_type::size_type size_type;
2383  typedef typename matrix_type::difference_type difference_type;
2384 
2385  for( difference_type c = A.cols( ) - 1 ; c >= 0 ; c-- )
2386  {
2387  double sum = b[ c ];
2388  for( size_type r = c + 1 ; r < A.rows( ) ; r++ )
2389  {
2390  sum -= A( r, c ) * b[ r ];
2391  }
2392 
2393  b[ c ] = sum / A( c, c );
2394  }
2395  }
2396 
2397 
2398  template < class Matrix >
2399  inline bool cholesky_factorization( const Matrix &H, Matrix &L, double lambda, double &lambda_modified )
2400  {
2401  typedef Matrix matrix_type;
2402  typedef typename matrix_type::value_type value_type;
2403  typedef typename matrix_type::size_type size_type;
2404  typedef typename matrix_type::difference_type difference_type;
2405 
2406  L.fill( 0 );
2407  for( size_type r = 0 ; r < H.rows( ) ; r++ )
2408  {
2409  double scale = H( r, r ) + lambda;
2410 
2411  for( size_type c = 0 ; c < r ; c++ )
2412  {
2413  scale -= L( r, c ) * L( r, c );
2414  }
2415 
2416  if( scale <= 0 )
2417  {
2418  // Rayleigh quotient trick を使ってλの補正値を計算する
2419  // CONDORの論文の 4.8 章
2420  matrix_type v( L.rows( ), 1 );
2421  v[ r ] = 1.0;
2422  for( difference_type j = r - 1 ; j >= 0 ; j-- )
2423  {
2424  double sum = 0.0;
2425  for( size_type i = j + 1 ; i <= r ; i++ )
2426  {
2427  sum += L( i, j ) * v[ i ];
2428  }
2429  v[ j ] = - sum / L( j, j );
2430  }
2431 
2432  lambda_modified = - scale / frobenius_norm( v );
2433 
2434  return( false );
2435  }
2436 
2437  scale = std::sqrt( scale );
2438  L( r, r ) = scale;
2439 
2440  for( size_type c = r + 1 ; c < H.cols( ) ; c++ )
2441  {
2442  double val = H( r, c );
2443  for( size_type l = 0 ; l < r ; l++ )
2444  {
2445  val -= L( c, l ) * L( r, l );
2446  }
2447 
2448  L( c, r ) = val / scale;
2449  }
2450  }
2451 
2452  return( true );
2453  }
2454 
2455  template < class Matrix >
2456  inline bool cholesky_factorization( const Matrix &H, Matrix &L, double lambda )
2457  {
2458  double dmy;
2459  return( cholesky_factorization( H, L, lambda, dmy ) );
2460  }
2461 
2462  template < class Matrix >
2463  inline bool compute_eigen_vector( const Matrix &L, Matrix &w, double lambda )
2464  {
2465  typedef Matrix matrix_type;
2466  typedef typename matrix_type::value_type value_type;
2467  typedef typename matrix_type::size_type size_type;
2468  typedef typename matrix_type::difference_type difference_type;
2469 
2470  w.resize( L.rows( ), 1 );
2471  for( size_type r = 0 ; r < L.rows( ) ; r++ )
2472  {
2473  if( L( r, r ) >= 0.0 )
2474  {
2475  w[ r ] = 1.0;
2476  }
2477  else
2478  {
2479  w[ r ] = -1.0;
2480  }
2481  }
2482 
2483  solve_( L, w );
2484  w *= 1.0 / frobenius_norm( w );
2485 
2486  return( true );
2487  }
2488 
2489  class polynomial : public matrix< double >
2490  {
2491  public:
2492  typedef matrix< double > base;
2493  typedef matrix< double > matrix_type;
2494  typedef matrix_type::value_type value_type;
2495  typedef matrix_type::size_type size_type;
2496  typedef matrix_type::difference_type difference_type;
2497  typedef matrix< difference_type > imatrix_type;
2498 
2499  private:
2500  size_type dimension;
2501  size_type N;
2502  imatrix_type alpha;
2503  imatrix_type alpha_;
2504  imatrix_type tr;
2505  imatrix_type tc;
2506  std::vector< double > r;
2507 
2508  public:
2509  polynomial( ) : dimension( 0 ), N( 1 )
2510  {
2511  }
2512 
2513  polynomial( size_type ndim ) : base( ( ndim + 1 ) * ( ndim + 2 ) / 2, 1 ), dimension( ndim ), N( ( ndim + 1 ) * ( ndim + 2 ) / 2 ), alpha( N, ndim ), alpha_( N, ndim ), tr( N, 1 ), tc( N, 1 ), r( ndim )
2514  {
2515  // 多項式補間を行うためのデータを生成する
2516  compute_polynomial_indeces( );
2517  }
2518 
2519  polynomial( const polynomial &poly ) : base( poly ), dimension( poly.dimension ), N( poly.N ), alpha( poly.alpha ), alpha_( poly.alpha_ ), tr( poly.tr ), tc( poly.tc ), r( poly.r )
2520  {
2521  }
2522 
2523  void reinitialize_polynomial( size_type ndim )
2524  {
2525  // 多項式を計算するためのデータをコピーする
2526  dimension = ndim;
2527  N = ( ndim + 1 ) * ( ndim + 2 ) / 2;
2528  alpha.resize( N, ndim );
2529  alpha_.resize( N, ndim );
2530  tr.resize( N, 1 );
2531  tc.resize( N, 1 );
2532  r.resize( ndim );
2533 
2534  base::resize( N, 1 );
2535 
2536  // 多項式補間を行うためのデータを生成する
2537  compute_polynomial_indeces( );
2538  }
2539 
2540  const polynomial & operator =( const polynomial &poly )
2541  {
2542  if( &poly == this )
2543  {
2544  return( *this );
2545  }
2546 
2547  // 基底クラスのデータをコピーする
2548  base::operator =( poly );
2549 
2550  // 多項式を計算するためのデータをコピーする
2551  dimension = poly.dimension;
2552  N = poly.N;
2553  alpha = poly.alpha;
2554  alpha_ = poly.alpha_;
2555  tr = poly.tr;
2556  tc = poly.tc;
2557  r = poly.r;
2558 
2559  return( *this );
2560  }
2561 
2562  const polynomial & operator =( const matrix_type &mat )
2563  {
2564  // 基底クラスのデータをコピーする
2565  base::operator =( mat );
2566 
2567  return( *this );
2568  }
2569 
2570  template < class Matrix >
2571  double operator ()( const Matrix &x )
2572  {
2573  difference_type n = dimension;
2574  const matrix_type &c = *this; // 多項式の係数ベクトル
2575 
2576  double r0 = c[ tr[ 0 ] ];
2577  for( difference_type j = 0 ; j < n ; j++ )
2578  {
2579  r[ j ] = 0.0;
2580  }
2581 
2582  for( size_type i = 1 ; i < N ; i++ )
2583  {
2584  difference_type k = tc[ i - 1 ];
2585  double rsum = r0;
2586  for( difference_type j = k ; j < n ; j++ )
2587  {
2588  rsum += r[ j ];
2589  r[ j ] = 0.0;
2590  }
2591 
2592  r0 = c[ tr[ i ] ];
2593  r[ k ] = x[ k ] * rsum;
2594  }
2595 
2596  double rsum = r0;
2597  for( difference_type j = 0 ; j < n ; j++ )
2598  {
2599  rsum += r[ j ];
2600  }
2601 
2602  return( rsum );
2603  }
2604 
2605  matrix_type compute_hessian_matrix( ) const
2606  {
2607  matrix_type H( dimension, dimension );
2608  const matrix_type &coeff = *this; // 多項式の係数ベクトル
2609 
2610  size_type indx = dimension + 1;
2611  for( size_type r = 0 ; r < H.rows( ) ; r++ )
2612  {
2613  H( r, r ) = coeff[ indx++ ] * 2.0;
2614  for( size_type c = r + 1 ; c < H.cols( ) ; c++ )
2615  {
2616  H( r, c ) = H( c, r ) = coeff[ indx++ ];
2617  }
2618  }
2619 
2620  return( H );
2621  }
2622 
2623  // 多項式の原点を平行移動する
2624  void translate( const matrix_type &x )
2625  {
2626  // 3.4.6 章参照
2627  polynomial &c = *this;
2628  matrix_type Hx = compute_hessian_matrix( ) * x;
2629 
2630  c[ 0 ] = c( x );
2631  for( size_type r = 0 ; r < Hx.rows( ) ; r++ )
2632  {
2633  c[ r + 1 ] += Hx[ r ];
2634  }
2635  }
2636 
2637 
2638  protected:
2639  void compute_polynomial_indeces( )
2640  {
2641  difference_type n = dimension;
2642  difference_type R = N;
2643  difference_type C = n;
2644  imatrix_type &lex = alpha_;
2645  imatrix_type &deg = alpha;
2646 
2647  lex.fill( 0 );
2648  deg.fill( 0 );
2649 
2650  {
2651  difference_type row = 1;
2652  for( difference_type c = 0 ; c < C ; c++ )
2653  {
2654  deg( row, c ) = 1;
2655  row++;
2656  }
2657 
2658  difference_type col = 0;
2659  difference_type val = 2;
2660  for( ; row < R && col < C ; )
2661  {
2662  switch( val )
2663  {
2664  case 2:
2665  deg( row, col ) = static_cast< imatrix_type::value_type >( val );
2666  val--;
2667  row++;
2668  break;
2669 
2670  case 1:
2671  for( difference_type c = 1 ; c < C - col ; c++ )
2672  {
2673  deg( row, col ) = static_cast< imatrix_type::value_type >( val );;
2674  deg( row, col + c ) = static_cast< imatrix_type::value_type >( val );;
2675  row++;
2676  }
2677  val--;
2678  col++;
2679  break;
2680 
2681  case 0:
2682  default:
2683  val = 2;
2684  break;
2685  }
2686  }
2687  }
2688 
2689  {
2690  difference_type col = 0;
2691  difference_type val = 2;
2692  for( difference_type row = 0 ; row < R - 1 ; )
2693  {
2694  switch( val )
2695  {
2696  case 2:
2697  lex( row, col ) = static_cast< imatrix_type::value_type >( val );
2698  val--;
2699  row++;
2700  break;
2701 
2702  case 1:
2703  for( difference_type c = 1 ; c < C - col ; c++ )
2704  {
2705  lex( row, col ) = static_cast< imatrix_type::value_type >( val );
2706  lex( row, col + c ) = static_cast< imatrix_type::value_type >( val );
2707  row++;
2708  }
2709  lex( row, col ) = static_cast< imatrix_type::value_type >( val );
2710  row++;
2711  val--;
2712  col++;
2713  break;
2714 
2715  case 0:
2716  default:
2717  val = 2;
2718  break;
2719  }
2720  }
2721 
2722  for( difference_type r = 0 ; r < R ; r++ )
2723  {
2724  difference_type c = C - 1;
2725  for( ; c >= 0 ; c-- )
2726  {
2727  if( lex( r, c ) != 0 )
2728  {
2729  break;
2730  }
2731  }
2732  tc[ r ] = static_cast< imatrix_type::value_type >( c );
2733  }
2734  }
2735 
2736  for( difference_type r = 0 ; r < R ; r++ )
2737  {
2738  difference_type indx = 0;
2739  for( ; indx < R ; indx++ )
2740  {
2741  difference_type c = 0;
2742  for( ; c < C ; c++ )
2743  {
2744  if( lex( r, c ) != deg( indx, c ) )
2745  {
2746  // 不一致
2747  break;
2748  }
2749  }
2750 
2751  if( c == C )
2752  {
2753  // 一致しているものを発見
2754  break;
2755  }
2756  }
2757 
2758  tr[ r ] = static_cast< imatrix_type::value_type >( indx );
2759  }
2760  }
2761  };
2762  }
2763 
2764 
2775  template < class Matrix, class Functor >
2776  size_t generate_first_point_set( const Matrix &xbase, std::vector< Matrix > &x, Matrix &f, Functor func, double rho )
2777  {
2778  typedef Matrix matrix_type;
2779  typedef typename matrix_type::value_type value_type;
2780  typedef typename matrix_type::size_type size_type;
2781  typedef typename matrix_type::difference_type difference_type;
2782 
2783  difference_type n = xbase.rows( );
2784  difference_type N = ( n + 1 ) * ( n + 2 ) / 2;
2785 
2786  // データの格納先を用意する
2787  f.resize( N, 1 );
2788  x.resize( N );
2789  for( difference_type i = 0 ; i < N ; i++ )
2790  {
2791  x[ i ].resize( n, 1 );
2792  }
2793 
2794  x[ 0 ] = xbase;
2795  f[ 0 ] = func( x[ 0 ] );
2796 
2797  for( difference_type j = 0 ; j < n ; j++ )
2798  {
2799  size_type k = j + 1;
2800  x[ k ] = xbase;
2801  x[ k ][ j ] += rho;
2802  f[ k ] = func( x[ k ] );
2803  }
2804 
2805  matrix_type s( n, 1 );
2806  for( difference_type j = 0 ; j < n ; j++ )
2807  {
2808  if( f[ j + 1 ] > f[ 0 ] )
2809  {
2810  s[ j ] = -1.0;
2811  }
2812  else
2813  {
2814  s[ j ] = +1.0;
2815  }
2816  }
2817 
2818  for( difference_type j = 0 ; j < n ; j++ )
2819  {
2820  difference_type k = j + 1 + n;
2821  if( s[ j ] < 0.0 )
2822  {
2823  x[ k ] = xbase;
2824  x[ k ][ j ] -= rho;
2825  }
2826  else
2827  {
2828  x[ k ] = xbase;
2829  x[ k ][ j ] += rho + rho;
2830  }
2831  f[ k ] = func( x[ k ] );
2832  }
2833 
2834  difference_type k = 2 * n + 1;
2835  for( difference_type j = 0 ; j < n ; j++ )
2836  {
2837  for( difference_type i = 0 ; i < j ; i++ )
2838  {
2839  x[ k ] = xbase;
2840  x[ k ][ i ] += rho * s[ i ];
2841  x[ k ][ j ] += rho * s[ j ];
2842  f[ k ] = func( x[ k ] );
2843  k++;
2844  }
2845  }
2846 
2847  // 評価関数値の小さい順に並べ替える
2848  for( size_type i = 0 ; i < x.size( ) ; i++ )
2849  {
2850  for( size_type j = 0 ; j < x.size( ) ; j++ )
2851  {
2852  if( f[ i ] < f[ j ] )
2853  {
2854  double tmp = f[ i ];
2855  f[ i ] = f[ j ];
2856  f[ j ] = tmp;
2857  x[ i ].swap( x[ j ] );
2858  }
2859  }
2860  }
2861 
2862  return( 0 );
2863  }
2864 
2874  template < class Matrix >
2875  bool compute_polynomial_basis( std::vector< Matrix > &x, Matrix &f, std::vector< __condor_utility__::polynomial > &poly_bases, __condor_utility__::polynomial &poly )
2876  {
2877  typedef Matrix matrix_type;
2878  typedef typename matrix_type::value_type value_type;
2879  typedef typename matrix_type::size_type size_type;
2880  typedef typename matrix_type::difference_type difference_type;
2881  typedef __condor_utility__::polynomial polynomial_type;
2882 
2883  if( x.size( ) != poly_bases.size( ) )
2884  {
2885  // 入力データと多項式の数が一致しない
2886  return( false );
2887  }
2888 
2889  size_type N = poly.size( );
2890  for( size_type k = 0 ; k < N ; k++ )
2891  {
2892  polynomial_type &pk = poly_bases[ k ];
2893 
2894  {
2895  // 多項式の正規化
2896  size_type index = k;
2897  double max = 0.0;
2898  if( k == 0 )
2899  {
2900  max = pk( x[ k ] );
2901  }
2902  else
2903  {
2904  for( size_type i = k ; i < N ; i++ )
2905  {
2906  double val = pk( x[ i ] );
2907  if( std::abs( val ) > std::abs( max ) )
2908  {
2909  max = val;
2910  index = i;
2911  }
2912 
2913  if( std::abs( max ) > 1.0 )
2914  {
2915  break;
2916  }
2917  }
2918  }
2919 
2920  // 分母が最大となるデータを用いる
2921  if( index != k )
2922  {
2923  // データを入れ替える
2924  x[ k ].swap( x[ index ] );
2925 
2926  double tmp = f[ k ];
2927  f[ k ] = f[ index ];
2928  f[ index ] = tmp;
2929  }
2930 
2931  pk /= max;
2932  }
2933 
2934  for( size_type j = 0 ; j < N ; j++ )
2935  {
2936  if( j != k )
2937  {
2938  polynomial_type &pj = poly_bases[ j ];
2939  pj -= pj( x[ k ] ) * pk;
2940  }
2941  }
2942  }
2943 
2944 
2945  // 最終的な多項式を求める
2946  poly.fill( 0 );
2947  for( size_type i = 0 ; i < poly_bases.size( ) ; i++ )
2948  {
2949  poly += f[ i ] * poly_bases[ i ];
2950  }
2951 
2952  return( true );
2953  }
2954 
2961  template < class Matrix >
2962  inline double compute_lambda_function1( const Matrix &H )
2963  {
2964  typedef Matrix matrix_type;
2965  typedef typename matrix_type::size_type size_type;
2966 
2967  double Hmin = H( 0, 0 );
2968  for( size_type i = 1 ; i < H.rows( ) ; i++ )
2969  {
2970  if( Hmin > H( i, i ) )
2971  {
2972  Hmin = H( i, i );
2973  }
2974  }
2975 
2976  return( -Hmin );
2977  }
2978 
2986  template < class Matrix >
2987  inline double compute_lambda_function2( const Matrix &H, double alpha = 1.0 )
2988  {
2989  typedef Matrix matrix_type;
2990  typedef typename matrix_type::size_type size_type;
2991 
2992  double max = -1.0e100;
2993  for( size_type i = 0 ; i < H.rows( ) ; i++ )
2994  {
2995  double val = H( i, i ) * alpha;
2996  for( size_type j = 0 ; j < i ; j++ )
2997  {
2998  val += std::abs( H( i, j ) );
2999  }
3000  for( size_type j = i + 1 ; j < H.cols( ) ; j++ )
3001  {
3002  val += std::abs( H( i, j ) );
3003  }
3004 
3005  if( max < val )
3006  {
3007  max = val;
3008  }
3009  }
3010 
3011  return( max );
3012  }
3013 
3014 
3022  template < class Matrix >
3023  inline double compute_lambda_lower_bound( const Matrix &H, double dg_delta )
3024  {
3025  double Hm = compute_lambda_function1( H );
3026  double Hg = compute_lambda_function2( H, +1.0 );
3027  double Hf = __condor_utility__::frobenius_norm( H );
3028  double Hi = __condor_utility__::infinitum_norm( H );
3029 
3030  return( __condor_utility__::maximum( 0, Hm, dg_delta - __condor_utility__::minimum( Hg, Hf, Hi ) ) );
3031  }
3032 
3033 
3041  template < class Matrix >
3042  inline double compute_lambda_upper_bound( const Matrix &H, double dg_delta )
3043  {
3044  double Hg = compute_lambda_function2( H, -1.0 );
3045  double Hf = __condor_utility__::frobenius_norm( H );
3046  double Hi = __condor_utility__::infinitum_norm( H );
3047 
3048  return( __condor_utility__::maximum( 0, dg_delta + __condor_utility__::minimum( Hg, Hf, Hi ) ) );
3049  }
3050 
3051 
3060  template < class Matrix >
3061  void compute_trust_region_step( const Matrix &xbest, Matrix &s, __condor_utility__::polynomial &poly, double delta, size_t max_loop = 1000 )
3062  {
3063  typedef Matrix matrix_type;
3064  typedef typename matrix_type::value_type value_type;
3065  typedef typename matrix_type::size_type size_type;
3066  typedef typename matrix_type::difference_type difference_type;
3067  typedef __condor_utility__::polynomial polynomial_type;
3068 
3069  matrix_type H = poly.compute_hessian_matrix( );
3070  matrix_type g = H * xbest;
3071  for( size_type i = 0 ; i < g.rows( ) ; i++ )
3072  {
3073  g[ i ] += poly[ i + 1 ];
3074  }
3075 
3076  double gnorm = __condor_utility__::frobenius_norm( g );
3077  double lambda = gnorm / delta;
3078  double lambdaL = compute_lambda_lower_bound( H, lambda );
3079  double lambdaU = compute_lambda_upper_bound( H, lambda );
3080 
3081  // lambdaL <= lambda <= lambdaU となるようにする
3082  lambda = __condor_utility__::maximum( lambdaL, __condor_utility__::minimum( lambda, lambdaU ) );
3083 
3084  bool cholesky_factorization_finished = false;
3085  matrix_type L( H.rows( ), H.cols( ) ), w( g.rows( ), g.cols( ) );
3086 
3087  // Trust Region step のベクトルを初期化する
3088  s.resize( g.rows( ), 1 );
3089  s.fill( 0 );
3090 
3091  size_type loop;
3092  for( loop = 0 ; loop < max_loop ; loop++ )
3093  {
3094  if( !cholesky_factorization_finished )
3095  {
3096  double lambdaM = 0.0;
3097  if( !__condor_utility__::cholesky_factorization( H, L, lambda, lambdaM ) )
3098  {
3099  lambdaL = __condor_utility__::maximum( lambdaL, lambda + lambdaM );
3100  lambda = __condor_utility__::maximum( std::sqrt( lambdaL * lambdaU ), lambdaL + 0.01 * ( lambdaU - lambdaL ) );
3101  continue;
3102  }
3103  }
3104  else
3105  {
3106  cholesky_factorization_finished = false;
3107  }
3108 
3109  s = -g;
3110 
3111  // L * L.t( ) * s = -g を解く
3112  // L が下三角行列であるということを利用して以下の手順で s を求める
3113  __condor_utility__::solve( L, s ); // まず L * ( L.t( ) * s' ) = -g を解く
3114  __condor_utility__::solve_( L, s ); // 次に L.t( ) * s = s' を解く
3115  double snorm = __condor_utility__::frobenius_norm( s );
3116 
3117  // 4.9.1 節の終了判定
3118  if( std::abs( snorm - delta ) < 0.01 * delta )
3119  {
3120  // 十分良い値が求まったと判定する
3121  s *= delta / snorm;
3122  break;
3123  }
3124  else if( snorm < delta )
3125  {
3126  if( lambda == 0.0 )
3127  {
3128  // 十分良い値が求まったと判定する
3129  break;
3130  }
3131 
3132  lambdaU = __condor_utility__::minimum( lambdaU, lambda );
3133 
3134 #if defined( __CHECK_HARD_CASE__ ) && __CHECK_HARD_CASE__ != 0
3135  // 複雑な場合かをチェックする
3136  matrix_type u;
3137  __condor_utility__::compute_eigen_vector( L, u, lambda );
3138 
3139  double uHu = __condor_utility__::inner_product( u, H, u, lambda );
3140 
3141  // Powell の UOBYQA アルゴリズムでの Hard Case 判定条件
3142  if( uHu <= 0.01 * ( uHu + lambda * delta * delta ) )
3143  {
3144  double a = 0.0;
3145  double b = 0.0;
3146  double c = - delta * delta;
3147  for( size_type i = 0 ; i < u.rows( ) ; i++ )
3148  {
3149  a += u[ i ] * u[ i ];
3150  b += s[ i ] * u[ i ];
3151  c += s[ i ] * s[ i ];
3152  }
3153 
3154  double alpha = 0.0;
3155  double bac = std::sqrt( b * b - a * c );
3156 
3157  double alpha1 = ( -b - bac ) / a;
3158  double alpha2 = ( -b + bac ) / a;
3159 
3160  matrix_type s1( s ), s2( s ), x1( xbest ), x2( xbest );
3161  for( size_type r = 0 ; r < s.rows( ) ; r++ )
3162  {
3163  s1[ r ] += alpha1 * u[ r ];
3164  s2[ r ] += alpha2 * u[ r ];
3165  x1[ r ] += s1[ r ];
3166  x2[ r ] += s2[ r ];
3167  }
3168 
3169  double f1 = poly( x1 );
3170  double f2 = poly( x2 );
3171 
3172  if( f1 >= f2 )
3173  {
3174  alpha = alpha1;
3175  s = s1;
3176  }
3177  else
3178  {
3179  alpha = alpha2;
3180  s = s2;
3181  }
3182 
3183  // Hard Case なので終了する
3184  break;
3185  }
3186 #endif
3187  }
3188  else// if( snorm > delta )
3189  {
3190  lambdaL = __condor_utility__::maximum( lambdaL, lambda );
3191  }
3192 
3193  w = s;
3194  __condor_utility__::solve( L, w );
3195 
3196  double wnorm = __condor_utility__::frobenius_norm( w );
3197  double lambdaN = lambdaL;
3198 
3199  // wnorm が限りなくゼロに近い場合への対処
3200  if( wnorm > 1.0e-20 )
3201  {
3202  lambdaN = lambda + ( ( snorm - delta ) / delta ) * ( ( snorm * snorm ) / ( wnorm * wnorm ) );
3203  lambdaN = __condor_utility__::maximum( lambdaL, __condor_utility__::minimum( lambdaN, lambdaU ) ); // 新しく求めたλの範囲をチェックする
3204  }
3205  else
3206  {
3207  lambdaN = __condor_utility__::maximum( lambdaL, lambdaU ); // 新しく求めたλの範囲をチェックする
3208  }
3209 
3210  // 新しく求めたλを使ってコレスキー分解を再度行う
3211  if( __condor_utility__::cholesky_factorization( H, L, lambdaN ) )
3212  {
3213  lambda = lambdaN;
3214  cholesky_factorization_finished = true;
3215  }
3216  else
3217  {
3218  lambdaL = __condor_utility__::maximum( lambdaL, lambdaN );
3219  lambda = __condor_utility__::maximum( std::sqrt( lambdaL * lambdaU ), lambdaL + 0.01 * ( lambdaU - lambdaL ) );
3220  }
3221  }
3222  }
3223 
3224 
3245  template < class T, class Allocator, class Functor >
3246  double minimization( matrix< T, Allocator > &p, Functor func, double rho, double rho_end, double tolerance, size_t &iterations, size_t max_iterations = 1000 )
3247  {
3248  typedef matrix< T, Allocator > matrix_type;
3249  typedef typename matrix_type::value_type value_type;
3250  typedef typename matrix_type::size_type size_type;
3251  typedef typename matrix_type::difference_type difference_type;
3252  typedef matrix< difference_type > imatrix_type;
3253  typedef __condor_utility__::polynomial polynomial_type;
3254 
3255  matrix_type f;
3256  std::vector< matrix_type > x;
3257 
3258  difference_type n = p.rows( );
3259  difference_type N = ( n + 1 ) * ( n + 2 ) / 2;
3260 
3261  // [Step 1] 2次のラグランジェ補間（超曲面）を計算するための初期点列を求める
3262  difference_type best_index = generate_first_point_set( p, x, f, func, rho );
3263 
3264 
3265  // [Step 2] ラグランジェ補間の基準となる点を決め，ラグランジェ多項式の係数を求める
3266  // 多項式を扱うオブジェクトを生成する
3267  // 最初のひとつを生成した後はすべてコピーする
3268  matrix_type xbase = x[ best_index ];
3269  polynomial_type poly;
3270  std::vector< polynomial_type > poly_bases( N );
3271  poly.reinitialize_polynomial( n );
3272 
3273  for( difference_type i = 0 ; i < N ; i++ )
3274  {
3275  poly_bases[ i ] = poly;
3276  }
3277  for( difference_type i = 0 ; i < N ; i++ )
3278  {
3279  poly_bases[ i ].fill( 0 );
3280  poly_bases[ i ][ i ] = 1.0;
3281 
3282  x[ i ] -= xbase;
3283  }
3284 
3285  // 多項式の基底を初期化する
3286  compute_polynomial_basis( x, f, poly_bases, poly );
3287 
3288 
3289  double M = 0.0; // 多項式補間のモデル評価値
3290  double delta = rho; // Trust Region の半径
3291  matrix_type tstep; // Trust Region 問題を解いた結果得られた関数値の減少方向
3292  double snorm = 0.0; // Trust Region 問題を解いた結果得られた関数値の減少方向の大きさ
3293  double Fold = f[ best_index ]; // 前回の反復操作で得られた関数値の最小値
3294  double Fold_ = Fold; // 許容相対誤差での終了判定用
3295  double Fnew = 1.0e100; // 新規に評価した関数値
3296  bool is_function_evaluated = false; // アルゴリズムの途中で，[Step 6]〜[Step 8]のスキップが発生したかどうか
3297  size_type number_of_updateM = 0;
3298 
3299  // 反復回数を初期化する
3300  iterations = 0;
3301 
3302  while( iterations < max_iterations )
3303  {
3304  Fold = __condor_utility__::minimum( Fnew, Fold );
3305  is_function_evaluated = false;
3306  for( size_type loop = 0 ; iterations++ < max_iterations ; loop++ )
3307  {
3308  double Fbest = f[ best_index ];
3309  const matrix_type &xbest = x[ best_index ];
3310 
3311  // [Step 5] Trust Region を求める
3312  compute_trust_region_step( xbest, tstep, poly, delta );
3313 
3314  // Trust Region から得られた新しい探索点を求める
3315  snorm = __condor_utility__::frobenius_norm( tstep );
3316  matrix_type xnew = xbest + tstep;
3317 
3318 
3319  // [Step 6] Trust Region ステップを終了し，モデルの再検証を行う
3320  // 1回は必ず以降の処理を実行するようにするため、初回はスキップする
3321  if( snorm < rho * 0.5 && loop > 0 )
3322  {
3323  break;
3324  }
3325 
3326 
3327  // [Step 7]
3328  // fbest を中心にラグランジェ多項式を当てはめるので，poly( xbest ) == fbest のはず
3329  double R = Fold - poly( xnew );
3330 
3331 
3332  // [Step 8] ノイズを付加して正しく最適化できているかを調べる？
3333  // 1回は必ず以降の処理を実行するようにするため、初回はスキップする
3334  double na = 0.0, nr = 0.0;
3335  double noise = 0.5 * __condor_utility__::maximum( na * ( 1.0 + nr ), nr * std::abs( Fbest ) );
3336  if( R < noise && loop > 0 )
3337  {
3338  // 今のところここにくることは無い
3339  break;
3340  }
3341 
3342 
3343  // [Step 9] 新しい位置で関数値を評価する
3344  {
3345  matrix_type ppp( xbase + xnew );
3346  Fnew = func( ppp );
3347  is_function_evaluated = true;
3348  }
3349 
3350 
3351  // [Step 10] モデルとの整合性を表す指標を計算する
3352  double r = R != 0.0 ? ( Fold - Fnew ) / R : Fold - Fnew;
3353 
3354 
3355  // [Step 11] Trust Region の半径を更新する
3356  if( 0.7 <= r )
3357  {
3358  delta = __condor_utility__::maximum( delta, snorm * 1.25, rho + snorm );
3359  }
3360  else if( 0.1 <= r )
3361  {
3362  delta = __condor_utility__::maximum( 0.5 * delta, snorm );
3363  }
3364  else
3365  {
3366  delta = 0.5 * snorm;
3367  }
3368 
3369  if( delta < 1.5 * rho )
3370  {
3371  delta = rho;
3372  }
3373 
3374 
3375  // [Step 12] 新しい点を補間点群に含め，多項式の係数を更新する
3376  // 3.4.3 と 3.4.4 を参照
3377  double model_step = 0.0;
3378  bool has_reduction = Fnew < Fbest;
3379  {
3380  difference_type index = -1;
3381  double max = -1.0e100;
3382  double Pnew = 0.0;
3383  if( Fnew < Fbest )
3384  {
3385  for( difference_type i = 0 ; i < N ; i++ )
3386  {
3387  double pnew = poly_bases[ i ]( xnew );
3388  double norm = __condor_utility__::frobenius_norm( x[ i ] - xnew ) / rho;
3389  double val = __condor_utility__::maximum( 1.0, norm * norm * norm ) * std::abs( pnew );
3390 
3391  if( val > max )
3392  {
3393  max = val;
3394  Pnew = pnew;
3395  index = i;
3396  }
3397  }
3398  }
3399  else
3400  {
3401  for( difference_type i = 0 ; i < N ; i++ )
3402  {
3403  if( i != best_index )
3404  {
3405  double pnew = poly_bases[ i ]( xnew );
3406  double norm = __condor_utility__::frobenius_norm( x[ i ] - xbest ) / rho;
3407  double val = __condor_utility__::maximum( 1.0, norm * norm * norm ) * std::abs( pnew );
3408 
3409  if( val > max )
3410  {
3411  max = val;
3412  Pnew = pnew;
3413  index = i;
3414  }
3415  }
3416  }
3417  }
3418 
3419  if( index >= 0 && Pnew != 0.0 )
3420  {
3421  // Model Step を計算する
3422  model_step = __condor_utility__::frobenius_norm( x[ index ] - xnew );
3423 
3424  // 多項式を更新する
3425  poly_bases[ index ] /= Pnew;
3426  for( difference_type i = 0 ; i < N ; i++ )
3427  {
3428  if( i != index )
3429  {
3430  poly_bases[ i ] -= poly_bases[ i ]( xnew ) * poly_bases[ index ];
3431  }
3432  }
3433 
3434 
3435  // データの置き換えを行う
3436  x[ index ] = xnew;
3437  f[ index ] = Fnew;
3438 
3439 
3440  // [Step 13] 最も良い値を持つものを更新する
3441  if( f[ best_index ] > Fnew )
3442  {
3443  Fold = __condor_utility__::minimum( Fnew, Fold );
3444  best_index = index;
3445  }
3446 
3447  // 最終的な多項式を求める
3448  poly.fill( 0 );
3449  for( size_type i = 0 ; i < poly_bases.size( ) ; i++ )
3450  {
3451  poly += f[ i ] * poly_bases[ i ];
3452  }
3453  }
3454  }
3455 
3456 
3457  // [Step 14] モデルの評価に用いる値 M を計算する
3458  {
3459  double sum = 0.0;
3460  for( size_type i = 0 ; i < poly_bases.size( ) ; i++ )
3461  {
3462  double len = __condor_utility__::frobenius_norm( xnew - x[ i ] );
3463  sum += std::abs( poly_bases[ i ]( xnew ) ) * len * len * len;
3464  }
3465 
3466  if( sum != 0.0 )
3467  {
3468  M = __condor_utility__::maximum( M, std::abs( poly( xnew ) - Fnew ) / sum );
3469  number_of_updateM++;
3470  }
3471  }
3472 
3473 
3474  // [Step 15] 関数の最小化に効果があったかどうかを判定する
3475  if( model_step > 2.0 * rho || snorm > 2.0 * rho || has_reduction )
3476  {
3477  }
3478  else
3479  {
3480  break;
3481  }
3482  }
3483 
3484  if( iterations >= max_iterations )
3485  {
3486  // 最大反復回数を経過した
3487  break;
3488  }
3489 
3490 
3491  // [Step 17.1] モデルチェックに使う評価基準値を求める
3492  double eps = 0.0;
3493  if( number_of_updateM < 10 )
3494  {
3495  eps = 0.0;
3496  }
3497  else if( snorm >= rho * 0.5 )
3498  {
3499  eps = 0.0;
3500  }
3501  else
3502  {
3503  // 4.10 章を参照
3504  // 多項式モデルの最も良い評価地の点付近の傾きを計算する
3505  matrix_type H = poly.compute_hessian_matrix( );
3506  matrix_type L( H.rows( ), H.cols( ) );
3507  double Hg = compute_lambda_function2( H, +1.0 );
3508  double Hf = __condor_utility__::frobenius_norm( H );
3509  double Hi = __condor_utility__::infinitum_norm( H );
3510 
3511  double lambdaL = 0.0;
3512  double lambdaU = __condor_utility__::minimum( Hg, Hf, Hi );
3513  while( lambdaL < 0.99 * lambdaU )
3514  {
3515  double lambda = ( lambdaL + lambdaU ) * 0.5;
3516 
3517  if( __condor_utility__::cholesky_factorization( H, L, -lambda ) )
3518  {
3519  lambdaL = lambda;
3520  }
3521  else
3522  {
3523  lambdaU = lambda;
3524  }
3525  }
3526 
3527  eps = 0.5 * rho * rho * ( lambdaL + lambdaU ) * 0.5;
3528  }
3529 
3530 
3531  // [17.2] モデルチェックに用いる点を求める（3.4.2章）
3532  {
3533  // まず，モデル評価に用いる点の集合Jを求める
3534  std::vector< __condor_utility__::__index_value_pair__ > J;
3535  J.reserve( N );
3536  for( difference_type i = 0 ; i < N ; i++ )
3537  {
3538  if( i != best_index )
3539  {
3540  double len = __condor_utility__::frobenius_norm( x[ i ] - x[ best_index ] );
3541  if( len > 2.0 * rho )
3542  {
3543  J.push_back( __condor_utility__::__index_value_pair__( i, len ) );
3544  }
3545  }
3546  }
3547 
3548  // モデルに悪影響を与える点の判定と多項式の更新を行う
3549  if( !J.empty( ) )
3550  {
3551  // 最も良い点からの距離値が大きい順に並べる
3552  std::sort( J.begin( ), J.end( ) );
3553 
3554  size_type index = 0;
3555  matrix_type d;
3556  for( ; index < J.size( ) ; index++ )
3557  {
3558  __condor_utility__::__index_value_pair__ &ivpair = J[ index ];
3559  polynomial_type &p = poly_bases[ ivpair.index ];
3560  matrix_type H = p.compute_hessian_matrix( );
3561  matrix_type g = H * x[ best_index ];
3562  for( size_type i = 0 ; i < g.rows( ) ; i++ )
3563  {
3564  g[ i ] += p[ i + 1 ];
3565  }
3566 
3567  double gnorm = __condor_utility__::frobenius_norm( g );
3568  if( gnorm == 0.0 )
3569  {
3570  continue;
3571  }
3572 
3573  matrix_type w( H.rows( ), 1 );
3574 
3575  {
3576  // ヘッセ行列 H の列ベクトルの乗る無が最大となるものを見つける
3577  // 式(5.5)と(5.6)を参照
3578  double max = -1.0e10;
3579  size_type col = 0;
3580  for( size_type c = 0 ; c < H.cols( ) ; c++ )
3581  {
3582  double sum = 0.0;
3583  for( size_type r = 0 ; r < H.rows( ) ; r++ )
3584  {
3585  sum += H( r, c ) * H( r, c );
3586  }
3587 
3588  if( max < sum )
3589  {
3590  max = sum;
3591  col = c;
3592  }
3593  }
3594 
3595  for( size_type r = 0 ; r < H.rows( ) ; r++ )
3596  {
3597  w[ r ] = H( r, col );
3598  }
3599  }
3600 
3601  matrix_type V = w;
3602  matrix_type D = H * w;
3603  double DHD = __condor_utility__::inner_product( D, H, D );
3604  double VD = __condor_utility__::inner_product( V, D );
3605  double DD = __condor_utility__::inner_product( D, D );
3606  double VHD = __condor_utility__::inner_product( V, H, D );
3607  double VHV = __condor_utility__::inner_product( V, H, V );
3608  double VV = __condor_utility__::inner_product( V, V );
3609 
3610  {
3611  double a = DHD * VD - DD * DD;
3612  double b = ( DHD * VV - DD * VD ) * 0.5;
3613  double c = DD * VV - VD * VD;
3614  double bac = std::sqrt( b * b - a * c );
3615  double r = 0.0;
3616 
3617  double r1 = ( -b - bac ) / a;
3618  double r2 = ( -b + bac ) / a;
3619 
3620  double f1 = ( r1 * r1 * DHD + 2.0 * r1 * VHD + VHV ) / ( VV + 2.0 * r1 * VD + r1 * r1 * DD );
3621  double f2 = ( r2 * r2 * DHD + 2.0 * r2 * VHD + VHV ) / ( VV + 2.0 * r2 * VD + r2 * r2 * DD );
3622 
3623  if( f1 > f2 )
3624  {
3625  r = r1;
3626  }
3627  else
3628  {
3629  r = r2;
3630  }
3631 
3632  // D を更新する
3633  D = V + r * D;
3634  }
3635 
3636  DD = __condor_utility__::inner_product( D, D );
3637  if( gnorm * DD + 2.0 * rho * std::abs( DHD ) < 0.5 )
3638  {
3639  double GD = __condor_utility__::inner_product( g, D );
3640  double scale = rho / std::sqrt( DD );
3641  if( GD * DHD < 0.0 )
3642  {
3643  d = - D * scale;
3644  }
3645  else
3646  {
3647  d = D * scale;
3648  }
3649  }
3650  else
3651  {
3652  // 5.2 章の u の組を作る
3653  // 大きさが1になるように正規化する
3654  // 式(5.10) からθを計算する
3655  matrix_type G = D / __condor_utility__::frobenius_norm( D );
3656  matrix_type V = g / gnorm;
3657  double VHG = __condor_utility__::inner_product( V, H, G );
3658  double VHV = __condor_utility__::inner_product( V, H, V );
3659  double GHG = __condor_utility__::inner_product( G, H, G );
3660  double theta = std::atan2( 2.0 * VHG, VHV - GHG ) * 0.5;
3661 
3662  // 式(5.9) を元に u を計算する
3663  double s_ = std::sin( theta );
3664  double c_ = std::cos( theta );
3665  matrix_type ut = c_ * G + s_ * V;
3666  matrix_type uh = -s_ * G + c_ * V;
3667 
3668  // 式(3.38)を最大にする d を得る
3669  // 5.3 章参照
3670  const double pai = 3.1415926535897932384626433832795;
3671  const double phi[ 4 ] = { 0.0, pai * 0.25, pai * 0.5, pai * 0.75 };
3672 
3673  // d の初期方向を設定する
3674  d = ( std::cos( phi[ 0 ] ) * uh + std::sin( phi[ 0 ] ) * ut );
3675  d = rho * d / __condor_utility__::frobenius_norm( d ); // 大きさを rho にする
3676 
3677  // 式(5.2)を最大にする d を選ぶ
3678  double max = std::abs( __condor_utility__::inner_product( g, d ) ) + 0.5 * std::abs( __condor_utility__::inner_product( d, H, d ) );
3679  for( size_type i = 1 ; i < 4 ; i++ )
3680  {
3681  matrix_type tmp = ( std::cos( phi[ i ] ) * uh + std::sin( phi[ i ] ) * ut );
3682  tmp = rho * tmp / __condor_utility__::frobenius_norm( tmp ); // 大きさを rho にする
3683  double val = std::abs( __condor_utility__::inner_product( g, tmp ) ) + 0.5 * std::abs( __condor_utility__::inner_product( tmp, H, tmp ) );
3684  if( val > max )
3685  {
3686  max = val;
3687  d = tmp;
3688  }
3689  }
3690 
3691  double GD = __condor_utility__::inner_product( g, d );
3692  double DHD = __condor_utility__::inner_product( d, H, d );
3693  if( GD * DHD < 0.0 )
3694  {
3695  d = - d;
3696  }
3697  }
3698 
3699  // 式(3.37)の評価
3700  double len = ivpair.value;
3701  double val = M * len * len * len * p( x[ best_index ] + d );
3702 
3703  if( val > eps )
3704  {
3705  break;
3706  }
3707  }
3708 
3709 
3710  // [17.3] 最もモデルに悪影響を与えると思われる点を入れ替える（3.4.4章）
3711  if( index < J.size( ) )
3712  {
3713  // 補間多項式の条件を満たさない点を発見
3714  difference_type outindex = J[ index ].index;
3715 
3716  matrix_type xnew = x[ best_index ] + d;
3717  matrix_type ppp = xbase + xnew;
3718  double fnew = func( ppp );
3719  is_function_evaluated = true;
3720 
3721  // 多項式を更新する
3722  poly_bases[ outindex ] /= poly_bases[ outindex ]( xnew );
3723  for( difference_type i = 0 ; i < N ; i++ )
3724  {
3725  if( i != outindex )
3726  {
3727  poly_bases[ i ] -= poly_bases[ i ]( xnew ) * poly_bases[ outindex ];
3728  }
3729  }
3730 
3731  // データの置き換えを行う
3732  x[ outindex ] = xnew;
3733  f[ outindex ] = fnew;
3734 
3735  // 最終的な多項式を求める
3736  poly.fill( 0 );
3737  for( size_type i = 0 ; i < poly_bases.size( ) ; i++ )
3738  {
3739  poly += f[ i ] * poly_bases[ i ];
3740  }
3741 
3742  // [Step 17.4] 最も良い値を持つものを更新する
3743  Fold = __condor_utility__::minimum( fnew, Fold );
3744  if( f[ best_index ] > fnew )
3745  {
3746  best_index = outindex;
3747  }
3748 
3749  // [Step 17.5] モデルが不適切と判断されたので値 M を再計算する
3750  {
3751  double sum = 0.0;
3752  for( size_type i = 0 ; i < poly_bases.size( ) ; i++ )
3753  {
3754  double len = __condor_utility__::frobenius_norm( xnew - x[ i ] );
3755  sum += std::abs( poly_bases[ i ]( xnew ) ) * len * len * len;
3756  }
3757 
3758  if( sum != 0.0 )
3759  {
3760  M = __condor_utility__::maximum( M, std::abs( poly( xnew ) - Fnew ) / sum );
3761  number_of_updateM++;
3762  }
3763  }
3764 
3765  // [Step 17.6] Step 4 へ戻る
3766  continue;
3767  }
3768  } // <- モデルに悪影響を与える点の判定と多項式の更新はここまで
3769 
3770  // [Step 17.7] モデルが正しく求まっている場合の判定
3771  if( snorm > rho )
3772  {
3773  continue;
3774  }
3775  }
3776 
3777  // [Step 18] ループの終了判定
3778  {
3779  if( rho <= rho_end || 2.0 * std::abs( Fold_ - Fold ) < tolerance * ( std::abs( Fold_ ) + std::abs( Fold ) ) )
3780  {
3781  break;
3782  }
3783  else
3784  {
3785  Fold_ = Fold;
3786  }
3787  }
3788 
3789  // [Step 19] Trust Region の半径を更新
3790  {
3791  double rho_old = rho;
3792  if( rho <= 16.0 * rho_end )
3793  {
3794  rho = rho_end;
3795  }
3796  else if( rho <= 250.0 * rho_end )
3797  {
3798  rho = std::sqrt( rho_end * rho );
3799  }
3800  else
3801  {
3802  rho *= 0.1;
3803  }
3804 
3805  delta = __condor_utility__::maximum( rho_old * 0.5, rho );
3806  }
3807 
3808 
3809  // [Step 20] 多項式の原点を平行移動する
3810  {
3811  matrix_type shift = x[ best_index ];
3812  xbase += x[ best_index ];
3813  poly.translate( shift );
3814  for( difference_type i = 0 ; i < N ; i++ )
3815  {
3816  poly_bases[ i ].translate( shift );
3817  if( i == best_index )
3818  {
3819  // 最も良い値を持つ者だけは 0 となるようにする
3820  x[ i ] *= 0.0;
3821  }
3822  else
3823  {
3824  x[ i ] -= shift;
3825  }
3826  }
3827  }
3828  }
3829 
3830 
3831  // [Step 21]
3832  {
3833  matrix_type xnew = xbase + x[ best_index ] + tstep;
3834 
3835  // 新しい位置で関数値を評価していない場合は評価を行う
3836  if( !is_function_evaluated )
3837  {
3838  Fnew = func( xnew );
3839  }
3840 
3841  if( Fold < Fnew )
3842  {
3843  p = xbase + x[ best_index ];
3844  }
3845  else
3846  {
3847  p = xnew;
3848  Fold = Fnew;
3849  }
3850  }
3851 
3852  return( Fold );
3853  }
3854 
3873  template < class T, class Allocator, class Functor >
3874  double minimization( matrix< T, Allocator > &p, Functor f, double rho, double rho_end, size_t max_iterations )
3875  {
3876  typedef __minimization_utility__::__no_copy_constructor_functor__< Functor > __no_copy_constructor_functor__;
3877  size_t itenum = 0;
3878  return( minimization( p, __no_copy_constructor_functor__( f ), rho, rho_end, 0.0, itenum, max_iterations ) );
3879  }
3880 
3899  template < class T, class Allocator, class Functor >
3900  double minimization( matrix< T, Allocator > &p, Functor f, double tolerance, size_t max_iterations )
3901  {
3902  typedef __minimization_utility__::__no_copy_constructor_functor__< Functor > __no_copy_constructor_functor__;
3903  size_t itenum = 0;
3904  return( minimization( p, __no_copy_constructor_functor__( f ), 1.0, 1.0e-8, tolerance, itenum, max_iterations ) );
3905  }
3906 }
3907 
3908 
3910 // 関数の最小化グループの終わり
3911 
3912 
3913 
3914 // mist名前空間の終わり
3915 _MIST_END
3916 
3917 
3918 #endif // __INCLUDE_MIST_MINIMIZATION__
3919 
Generated on Wed Nov 12 2014 19:44:19 for MIST by doxygen 1.8.1.2