- メインページ
- モジュール
- ネームスペース
- データ構造
- ファイル

- ファイル一覧
- グローバル

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