1 /*
2 * This source code is part of
3 *
4 * E R K A L E
5 * -
6 * HF/DFT from Hel
7 *
8 * Written by Susi Lehtola, 2010-2013
9 * Copyright (c) 2010-2013, Susi Lehtola
10 *
11 * This program is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU General Public License
13 * as published by the Free Software Foundation; either version 2
14 * of the License, or (at your option) any later version.
15 */
16
17 #include "unitary.h"
18 #include "timer.h"
19 #include "mathf.h"
20 #include "linalg.h"
21 #include <cfloat>
22
23 #ifdef SVNRELEASE
24 #include "version.h"
25 #endif
26
27 #ifdef _OPENMP
28 #include <omp.h>
29 #endif
30
UnitaryFunction(int qv,bool max)31 UnitaryFunction::UnitaryFunction(int qv, bool max): W(arma::cx_mat()), f(0.0), q(qv) {
32 /// Maximize or minimize?
33 sign = max ? 1 : -1;
34 }
35
~UnitaryFunction()36 UnitaryFunction::~UnitaryFunction() {
37 }
38
setW(const arma::cx_mat & Wv)39 void UnitaryFunction::setW(const arma::cx_mat & Wv) {
40 W=Wv;
41 }
42
getW() const43 arma::cx_mat UnitaryFunction::getW() const {
44 return W;
45 }
46
getq() const47 int UnitaryFunction::getq() const {
48 return q;
49 }
50
getf() const51 double UnitaryFunction::getf() const {
52 return f;
53 }
54
getsign() const55 int UnitaryFunction::getsign() const {
56 return sign;
57 }
58
converged()59 bool UnitaryFunction::converged() {
60 /// Dummy default function
61 return false;
62 }
63
legend() const64 std::string UnitaryFunction::legend() const {
65 /// Dummy default function
66 return "";
67 }
68
status(bool lfmt)69 std::string UnitaryFunction::status(bool lfmt) {
70 /// Dummy default function
71 (void) lfmt;
72 return "";
73 }
74
UnitaryOptimizer(double Gthrv,double Fthrv,bool ver,bool realv)75 UnitaryOptimizer::UnitaryOptimizer(double Gthrv, double Fthrv, bool ver, bool realv) : G(arma::cx_mat()), H(arma::cx_mat()), Hvec(arma::cx_mat()), Hval(arma::vec()), Tmu(0.0), verbose(ver), real(realv), Gthr(Gthrv), Fthr(Fthrv) {
76
77 // Defaults
78 // use 3rd degree polynomial to fit derivative
79 polynomial_degree=4;
80 // and in the fourier transform use
81 fourier_periods=5; // five quasi-periods
82 fourier_samples=3; // three points per period
83
84 debug=false;
85
86 // Logfile is closed
87 log=NULL;
88 }
89
~UnitaryOptimizer()90 UnitaryOptimizer::~UnitaryOptimizer() {
91 if(log!=NULL)
92 fclose(log);
93 }
94
set_debug(bool d)95 void UnitaryOptimizer::set_debug(bool d) {
96 debug=d;
97 }
98
set_thr(double Geps,double Feps)99 void UnitaryOptimizer::set_thr(double Geps, double Feps) {
100 Gthr=Geps;
101 Fthr=Feps;
102 }
103
open_log(const std::string & fname)104 void UnitaryOptimizer::open_log(const std::string & fname) {
105 if(log!=NULL)
106 fclose(log);
107
108 if(fname.length()) {
109 log=fopen(fname.c_str(),"w");
110 #ifdef _OPENMP
111 fprintf(log,"ERKALE - Localization from Hel, OpenMP version, running on %i cores.\n",omp_get_max_threads());
112 #else
113 fprintf(log,"ERKALE - Localization from Hel, serial version.\n");
114 #endif
115 fprint_copyright(log);
116 fprint_license(log);
117 #ifdef SVNRELEASE
118 fprintf(log,"At svn revision %s.\n\n",SVNREVISION);
119 #endif
120 fprint_hostname(log);
121 }
122 }
123
get_rotation(double step) const124 arma::cx_mat UnitaryOptimizer::get_rotation(double step) const {
125 // Rotation matrix is
126 arma::cx_mat rot=Hvec*arma::diagmat(arma::exp(step*COMPLEXI*Hval))*arma::trans(Hvec);
127 if(real)
128 // Zero out possible imaginary part
129 rot=arma::real(rot)*COMPLEX1;
130
131 return rot;
132 }
133
set_poly(int deg)134 void UnitaryOptimizer::set_poly(int deg) {
135 polynomial_degree=deg;
136 }
137
set_fourier(int samples,int pers)138 void UnitaryOptimizer::set_fourier(int samples, int pers) {
139 fourier_periods=pers;
140 fourier_samples=samples;
141 }
142
update_gradient(const arma::cx_mat & W,UnitaryFunction * f)143 void UnitaryOptimizer::update_gradient(const arma::cx_mat & W, UnitaryFunction *f) {
144 // Euclidean gradient
145 arma::cx_mat Gammak;
146 double J;
147 f->cost_func_der(W,J,Gammak);
148 // Fix sign
149 // Gammak*=f->getsign();
150
151 // Riemannian gradient, Abrudan 2009 table 3 step 2
152 G=Gammak*arma::trans(W) - W*arma::trans(Gammak);
153 }
154
update_search_direction(int q)155 void UnitaryOptimizer::update_search_direction(int q) {
156 // Diagonalize -iH to find eigenvalues purely imaginary
157 // eigenvalues iw_i of H; Abrudan 2009 table 3 step 1.
158 bool diagok=arma::eig_sym(Hval,Hvec,-COMPLEXI*H);
159 if(!diagok) {
160 ERROR_INFO();
161 throw std::runtime_error("Unitary optimization: error diagonalizing H.\n");
162 }
163
164 // Max step length
165 double wmax=arma::max(arma::abs(Hval));
166 Tmu=2.0*M_PI/(q*wmax);
167 }
168
optimize(UnitaryFunction * & f,enum unitmethod met,enum unitacc acc,size_t maxiter)169 double UnitaryOptimizer::optimize(UnitaryFunction* & f, enum unitmethod met, enum unitacc acc, size_t maxiter) {
170 // Get the matrix
171 arma::cx_mat W=f->getW();
172 if(real)
173 W=arma::real(W)*COMPLEX1;
174
175 // Current and old gradient
176 arma::cx_mat oldG;
177 G.zeros(W.n_cols,W.n_cols);
178 // Current and old search direction
179 arma::cx_mat oldH;
180 H.zeros(W.n_cols,W.n_cols);
181
182 if(W.n_cols<2) {
183 // No optimization is necessary.
184 W.eye();
185 return f->cost_func(W);
186 }
187
188 // Check matrix
189 ::check_unitarity(W);
190 // Check derivative
191 check_derivative(f);
192
193 // Info from previous iteration
194 UnitaryFunction *oldf=NULL;
195
196 // Iteration number
197 size_t k=0;
198
199 // Print out the legend
200 if(verbose)
201 print_legend(f);
202
203 if(maxiter>0)
204 while(true) {
205 Timer t;
206
207 // Store old gradient and search direction
208 oldG=G;
209 oldH=H;
210
211 // Compute the cost function and the euclidean derivative, Abrudan 2009 table 3 step 2
212 update_gradient(W,f);
213
214 // Compute update coefficient
215 double gamma=0.0;
216 if(acc==SDSA || (k-1)%W.n_cols==0) {
217 // Reset step in CG, or steepest descent / steepest ascent
218 gamma=0.0;
219 } else if(acc==CGPR) {
220 // Compute Polak-Ribière coefficient
221 gamma=bracket(G - oldG, G) / bracket(oldG, oldG);
222 } else if(acc==CGFR) {
223 // Fletcher-Reeves
224 gamma=bracket(G, G) / bracket(oldG, oldG);
225 } else if(acc==CGHS) {
226 // Hestenes-Stiefel
227 gamma=bracket(G - oldG, G) / bracket(G - oldG, oldH);
228 } else
229 throw std::runtime_error("Unsupported update.\n");
230
231 // Perform update
232 if(gamma==0.0) {
233 // Use gradient direction
234 H=G;
235 } else {
236 // Compute update
237 H=G+gamma*H;
238 // Make sure H stays skew symmetric
239 H=0.5*(H-arma::trans(H));
240
241 // Check that update is OK
242 if(bracket(G,H)<0.0) {
243 H=G;
244 printf("CG search direction reset.\n");
245 }
246 }
247 // Update search direction
248 update_search_direction(f->getq());
249
250 // Save old iteration data
251 if(oldf)
252 delete oldf;
253 oldf=f->copy();
254
255 // Take a step
256 if(met==POLY_DF) {
257 polynomial_step_df(f);
258 } else if(met==POLY_F) {
259 polynomial_step_f(f);
260 } else if(met==FOURIER_DF) {
261 fourier_step_df(f);
262 } else if(met==ARMIJO) {
263 armijo_step(f);
264 } else {
265 ERROR_INFO();
266 throw std::runtime_error("Method not implemented.\n");
267 }
268
269 // Increase iteration number
270 k++;
271 // Update matrix
272 W=f->getW();
273
274 // Print progress
275 if(verbose) {
276 print_progress(k,f,oldf);
277 print_time(t);
278 }
279
280 // Check for convergence. Don't do the check in the first
281 // iteration, because for canonical orbitals it can be a really
282 // bad saddle point
283 double J=f->getf();
284
285 double oldJ = (oldf==NULL) ? 0.0 : oldf->getf();
286 if( k>1 && (f->converged() || (bracket(G,G)<Gthr && fabs(J-oldJ)<Fthr))) {
287 if(verbose) {
288 printf("Converged.\n");
289 fflush(stdout);
290
291 // Print classification
292 classify(W);
293 }
294
295 break;
296 } else if(k==maxiter) {
297 if(verbose) {
298 printf(" %s\nNot converged.\n",t.elapsed().c_str());
299 fflush(stdout);
300 }
301
302 break;
303 }
304
305 if(debug) {
306 char fname[80];
307 sprintf(fname,"unitary_%04i.dat",(int) k);
308 FILE *p=fopen(fname,"w");
309 UnitaryFunction *uf=f->copy();
310 for(int i=-80;i<=80;i++) {
311 double x=i*0.05;
312 double xT=x*Tmu;
313
314 double y=uf->cost_func(get_rotation(xT)*W);
315 fprintf(p,"% e % e % e\n",x,xT,y);
316 }
317 fclose(p);
318 delete uf;
319 }
320 }
321
322 if(oldf)
323 delete oldf;
324
325 return f->getf();
326 }
327
print_legend(const UnitaryFunction * f) const328 void UnitaryOptimizer::print_legend(const UnitaryFunction *f) const {
329 printf(" %4s %13s %13s %12s %s\n","iter","J","delta J","<G,G>",f->legend().c_str());
330 }
331
print_progress(size_t k,UnitaryFunction * f,const UnitaryFunction * fold) const332 void UnitaryOptimizer::print_progress(size_t k, UnitaryFunction *f, const UnitaryFunction *fold) const {
333 double J=f->getf();
334 if(fold==NULL)
335 // No info on delta J
336 printf(" %4i % e %13s %e %s",(int) k,J,"",bracket(G,G),f->status().c_str());
337 else {
338 double oldJ=fold->getf();
339 printf(" %4i % e % e %e %s",(int) k,J,J-oldJ,bracket(G,G),f->status().c_str());
340 }
341 fflush(stdout);
342
343 if(log!=NULL) {
344 fprintf(log,"%4i % .16e %.16e %s",(int) k,J,bracket(G,G),f->status(true).c_str());
345 fflush(log);
346 }
347 }
348
print_time(const Timer & t) const349 void UnitaryOptimizer::print_time(const Timer & t) const {
350 printf(" %s\n",t.elapsed().c_str());
351 fflush(stdout);
352
353 if(log!=NULL) {
354 fprintf(log,"%e\n",t.get());
355 fflush(log);
356 }
357 }
358
print_step(enum unitmethod & met,double step) const359 void UnitaryOptimizer::print_step(enum unitmethod & met, double step) const {
360 /*
361 if(met==POLY_DF)
362 printf("Polynomial_df step %e (%e of Tmu)\n",step,step/Tmu);
363 else if(met==POLY_FDF)
364 printf("Polynomial_fdf step %e (%e of Tmu)\n",step,step/Tmu);
365 else if(met==FOURIER_DF)
366 printf("Fourier_df step %e (%e of Tmu)\n",step,step/Tmu);
367 else if(met==ARMIJO)
368 printf("Armijo step %e (%e of Tmu)\n",step,step/Tmu);
369 else {
370 ERROR_INFO();
371 throw std::runtime_error("Method not implemented.\n");
372 }
373 */
374
375 (void) met;
376 (void) step;
377 if(log!=NULL)
378 fprintf(log,"%e\n",step);
379 }
380
classify(const arma::cx_mat & W) const381 void UnitaryOptimizer::classify(const arma::cx_mat & W) const {
382 if(real)
383 return;
384
385 // Classify matrix
386 double realpart=rms_norm(arma::real(W));
387 double imagpart=rms_norm(arma::imag(W));
388
389 printf("Transformation matrix is");
390 if(imagpart<sqrt(DBL_EPSILON)*realpart)
391 printf(" real");
392 else if(realpart<sqrt(DBL_EPSILON)*imagpart)
393 printf(" imaginary");
394 else
395 printf(" complex");
396
397 printf(", re norm %e, im norm %e\n",realpart,imagpart);
398 }
399
check_derivative(const UnitaryFunction * fp0)400 void UnitaryOptimizer::check_derivative(const UnitaryFunction *fp0) {
401 UnitaryFunction *fp=fp0->copy();
402
403 arma::cx_mat W0=fp0->getW();
404
405 // Compute gradient
406 update_gradient(W0,fp);
407 // and the search direction
408 arma::cx_mat Hs=G;
409 update_search_direction(fp->getq());
410
411 // Get cost function and derivative matrix
412 arma::cx_mat der;
413 double Jo;
414 fp->cost_func_der(W0,Jo,der);
415 // The derivative is
416 double dfdmu=step_der(W0,der);
417
418 // Compute trial value.
419 double trstep=Tmu*sqrt(DBL_EPSILON);
420 arma::cx_mat Wtr=get_rotation(trstep*fp0->getsign())*W0;
421 double Jtr=fp->cost_func(Wtr);
422
423 // Estimated change in function is
424 double dfest=trstep*dfdmu;
425 // Real change in function is
426 double dfreal=Jtr-Jo;
427
428 // Is the difference ok? Check absolute or relative magnitude
429 if(fabs(dfest)>sqrt(DBL_EPSILON)*std::max(1.0,fabs(fp->getf())) && fabs(dfest-dfreal)>1e-2*fabs(dfest)) {
430 fprintf(stderr,"\nDerivative mismatch error!\n");
431 fprintf(stderr,"Used step size %e, value of function % e.\n",trstep,Jo);
432 fprintf(stderr,"Estimated change of function % e\n",dfest);
433 fprintf(stderr,"Realized change of function % e\n",dfreal);
434 fprintf(stderr,"Difference in changes % e\n",dfest-dfreal);
435 fprintf(stderr,"Relative error in changes % e\n",(dfest-dfreal)/dfest);
436 throw std::runtime_error("Derivative mismatch! Check your cost function and its derivative.\n");
437 }
438
439 delete fp;
440 }
441
polynomial_step_f(UnitaryFunction * & fp)442 void UnitaryOptimizer::polynomial_step_f(UnitaryFunction* & fp) {
443 // Amount of points to use is
444 int npoints=polynomial_degree;
445 // Spacing
446 double deltaTmu=Tmu/(npoints-1);
447 // Step size
448 double step=0.0;
449 arma::cx_mat W=fp->getW();
450
451 UnitaryFunction* fline[npoints];
452 for(int i=0;i<npoints;i++)
453 fline[i]=fp->copy();
454
455 while(step==0.0) {
456 // Evaluate the cost function at the expansion points
457 arma::vec mu(npoints);
458 arma::vec f(npoints);
459
460 // 0th point is current point!
461 for(int i=0;i<npoints;i++) {
462 // Mu in the point is
463 mu(i)=i*deltaTmu;
464
465 // Trial matrix is
466 arma::cx_mat Wtr=get_rotation(mu(i)*fp->getsign())*W;
467 // and the function is
468 f(i)=fline[i]->cost_func(Wtr);
469 }
470
471 // Fit to polynomial of order p
472 arma::vec coeff=fit_polynomial(mu,f);
473
474 // Find out zeros of the derivative polynomial
475 arma::vec roots=solve_roots(derivative_coefficients(coeff));
476 // get the smallest positive one
477 step=smallest_positive(roots);
478 if(step==0.0) {
479 deltaTmu/=2.0;
480 printf("No root found, halving step size to %e.\n",deltaTmu);
481 continue;
482 }
483
484 // Is the step length in the allowed region?
485 if(step>0.0 && step <=Tmu) {
486 // Yes. Calculate the new value
487 arma::cx_mat Wtr=get_rotation(step*fp->getsign())*W;
488 UnitaryFunction *newf=fp->copy();
489
490 double J=fp->getf();
491 double Jtr=newf->cost_func(Wtr);
492
493 // Accept the step?
494 if( fp->getsign()*(Jtr-J) > 0.0) {
495 // Yes.
496 delete fp;
497 fp=newf;
498 // Free memory
499 for(int i=0;i<npoints;i++)
500 delete fline[i];
501
502 return;
503 } else
504 delete newf;
505 }
506
507 // If we are still here, then just get the minimum value
508 if(step==0.0 || step > Tmu) {
509 fprintf(stderr,"Line search interpolation failed.\n");
510 fflush(stderr);
511 double minval=arma::max(fp->getsign()*f);
512 for(size_t i=0;i<mu.n_elem;i++)
513 if(minval==f(i)) {
514 delete fp;
515 fp=fline[i];
516 break;
517 }
518 }
519 }
520
521 for(int i=0;i<npoints;i++)
522 delete fline[i];
523
524 return;
525 }
526
polynomial_step_df(UnitaryFunction * & fp)527 void UnitaryOptimizer::polynomial_step_df(UnitaryFunction* & fp) {
528 // Matrix
529 arma::cx_mat W=fp->getW();
530 // Amount of points to use is
531 int npoints=polynomial_degree;
532 // Spacing
533 double deltaTmu=Tmu/(npoints-1);
534 int halved=0;
535
536 UnitaryFunction* fline[npoints];
537 for(int i=0;i<npoints;i++)
538 fline[i]=fp->copy();
539
540 while(true) {
541 // Evaluate the cost function at the expansion points
542 arma::vec mu(npoints);
543 arma::vec fd(npoints);
544 arma::vec fv(npoints);
545
546 for(int i=0;i<npoints;i++) {
547 // Mu in the point is
548 mu(i)=i*deltaTmu;
549
550 // Trial matrix is
551 arma::cx_mat Wtr=get_rotation(mu(i)*fp->getsign())*W;
552 // and the function is
553 arma::cx_mat der;
554 //der=fline[i]->cost_der(Wtr);
555 fline[i]->cost_func_der(Wtr,fv(i),der);
556 // and the derivative is
557 fd(i)=step_der(Wtr,der);
558 }
559
560 // Sanity check - is derivative of the right sign?
561 if(fd(0)<0.0) {
562 printf("Derivative is of the wrong sign!\n");
563 arma::trans(mu).print("mu");
564 arma::trans(fd).print("J'(mu)");
565 // throw std::runtime_error("Derivative consistency error.\n");
566 fprintf(stderr,"Warning - inconsistent sign of derivative.\n");
567 }
568
569 // Fit to polynomial of order p
570 arma::vec coeff=fit_polynomial(mu,fd);
571
572 // Find out zeros of the polynomial
573 arma::vec roots=solve_roots(coeff);
574 // get the smallest positive one
575 double step=smallest_positive(roots);
576
577 /*
578 printf("Trial step size is %e.\n",step);
579 mu.t().print("mu");
580 fd.t().print("f'");
581 fv.t().print("f");
582 */
583
584 // Is the step length in the allowed region?
585 if(step>0.0 && step <=Tmu) {
586 // Yes. Calculate the new value
587 arma::cx_mat Wtr=get_rotation(step*fp->getsign())*W;
588 UnitaryFunction *newf=fp->copy();
589
590 double J=fp->getf();
591 double Jtr=newf->cost_func(Wtr);
592
593 // Accept the step?
594 if( fp->getsign()*(Jtr-J) > 0.0) {
595 // Yes.
596 // printf("Function value changed by %e, accept.\n",Jtr-J);
597 delete fp;
598 fp=newf;
599 break;
600 } else {
601 fprintf(stderr,"Line search interpolation failed.\n");
602 fflush(stderr);
603 if(halved<1) {
604 delete newf;
605 halved++;
606 deltaTmu/=2.0;
607 continue;
608 } else {
609 // Try an Armijo step
610 return armijo_step(fp);
611
612 /*
613 ERROR_INFO();
614 throw std::runtime_error("Problem in polynomial line search - could not find suitable extremum!\n");
615 */
616 }
617 }
618
619 delete newf;
620 } else {
621 fprintf(stderr,"Line search interpolation failed.\n");
622 fflush(stderr);
623 if(halved<4) {
624 halved++;
625 deltaTmu/=2.0;
626 continue;
627 } else {
628 ERROR_INFO();
629 throw std::runtime_error("Problem in polynomial line search - could not find suitable extremum!\n");
630 }
631 }
632 }
633
634 for(int i=0;i<npoints;i++)
635 delete fline[i];
636 }
637
step_der(const arma::cx_mat & W,const arma::cx_mat & der) const638 double UnitaryOptimizer::step_der(const arma::cx_mat & W, const arma::cx_mat & der) const {
639 return 2.0*std::real(arma::trace(der*arma::trans(W)*arma::trans(H)));
640 }
641
armijo_step(UnitaryFunction * & fp)642 void UnitaryOptimizer::armijo_step(UnitaryFunction* & fp) {
643 // Start with half of maximum.
644 double step=Tmu/2.0;
645
646 // Initial rotation matrix
647 arma::cx_mat R=get_rotation(step*fp->getsign());
648
649 // Helper
650 UnitaryFunction *hlp=fp->copy();
651
652 // Original rotation
653 arma::cx_mat W(fp->getW());
654
655 // Current value
656 double J=fp->getf();
657
658 // Evaluate function at R2
659 double J2=hlp->cost_func(R*R*W);
660
661 if(fp->getsign()==-1) {
662 // Minimization.
663
664 // First condition: f(W) - f(R^2 W) >= mu*<G,H>
665 while(J-J2 >= step*bracket(G,H)) {
666 // Increase step size.
667 step*=2.0;
668 R=get_rotation(step*fp->getsign());
669
670 // and re-evaluate J2
671 J2=hlp->cost_func(R*R*W);
672 }
673
674 // Evaluate function at R
675 double J1=hlp->cost_func(R*W);
676
677 // Second condition: f(W) - f(R W) <= mu/2*<G,H>
678 while(J-J1 < step/2.0*bracket(G,H)) {
679 // Decrease step size.
680 step/=2.0;
681 R=get_rotation(step*fp->getsign());
682
683 // and re-evaluate J1
684 J1=hlp->cost_func(R*W);
685 }
686
687 } else if(fp->getsign()==1) {
688 // Maximization
689
690 // First condition: f(W) - f(R^2 W) >= mu*<G,H>
691 while(J-J2 <= -step*bracket(G,H)) {
692 // Increase step size.
693 step*=2.0;
694 R=get_rotation(step*fp->getsign());
695
696 // and re-evaluate J2
697 J2=hlp->cost_func(R*R*W);
698 }
699
700 // Evaluate function at R
701 double J1=hlp->cost_func(R*W);
702
703 // Second condition: f(W) - f(R W) <= mu/2*<G,H>
704 while(J-J1 > -step/2.0*bracket(G,H)) {
705 // Decrease step size.
706 step/=2.0;
707 R=get_rotation(step*fp->getsign());
708
709 // and re-evaluate J1
710 J1=hlp->cost_func(R*W);
711 }
712 } else
713 throw std::runtime_error("Invalid optimization direction!\n");
714
715 // Update solution
716 delete fp;
717 fp=hlp;
718 }
719
fourier_shift(const arma::cx_vec & c)720 arma::cx_vec fourier_shift(const arma::cx_vec & c) {
721 // Amount of elements
722 size_t N=c.n_elem;
723
724 // Midpoint is at at
725 size_t m=N/2;
726 if(N%2==1)
727 m++;
728
729 // Returned vector
730 arma::cx_vec ret(N);
731 ret.zeros();
732
733 // Low frequencies
734 ret.subvec(0,N-1-m)=c.subvec(m,N-1);
735 // High frequencies
736 ret.subvec(N-m,N-1)=c.subvec(0,m-1);
737
738 return ret;
739 }
740
fourier_step_df(UnitaryFunction * & f)741 void UnitaryOptimizer::fourier_step_df(UnitaryFunction* & f) {
742 // Length of DFT interval
743 double fourier_interval=fourier_periods*Tmu;
744 // and of the transform. We want integer division here!
745 int fourier_length=2*((fourier_samples*fourier_periods)/2)+1;
746
747 // Step length is
748 double deltaTmu=fourier_interval/fourier_length;
749
750 // Helpers
751 UnitaryFunction * fs[fourier_length];
752 for(int i=0;i<fourier_length;i++)
753 fs[i]=f->copy();
754 arma::cx_mat W=f->getW();
755
756 // Values of mu, J(mu) and J'(mu)
757 arma::vec mu(fourier_length);
758 arma::vec fv(fourier_length);
759 arma::vec fp(fourier_length);
760 for(int i=0;i<fourier_length;i++) {
761 // Value of mu is
762 mu(i)=i*deltaTmu;
763
764 // Trial matrix is
765 arma::cx_mat Wtr=get_rotation(mu(i)*f->getsign())*W;
766 arma::cx_mat der;
767 fs[i]->cost_func_der(Wtr,fv(i),der);
768
769 // Compute the derivative
770 fp(i)=step_der(Wtr,der);
771 }
772
773 // Compute Hann window
774 arma::vec hannw(fourier_length);
775 for(int i=0;i<fourier_length;i++)
776 hannw(i)=0.5*(1-cos((i+1)*2.0*M_PI/(fourier_length+1.0)));
777
778 // Windowed derivative is
779 arma::vec windowed(fourier_length);
780 for(int i=0;i<fourier_length;i++)
781 windowed(i)=fp(i)*hannw(i);
782
783 // Fourier coefficients
784 arma::cx_vec coeffs=arma::fft(windowed)/fourier_length;
785
786 // Reorder coefficients
787 arma::cx_vec shiftc=fourier_shift(coeffs);
788
789 // Find roots of polynomial
790 arma::cx_vec croots=solve_roots_cplx(shiftc);
791
792 // Figure out roots on the unit circle
793 double circletol=1e-2;
794 std::vector<double> muval;
795 for(size_t i=0;i<croots.n_elem;i++)
796 if(fabs(std::abs(croots(i))-1)<circletol) {
797 // Root is on the unit circle. Angle is
798 double phi=std::imag(std::log(croots(i)));
799
800 // Convert to the real length scale
801 phi*=fourier_interval/(2*M_PI);
802
803 // Avoid aliases
804 phi=fmod(phi,fourier_interval);
805 // and check for negative values (fmod can return negative values)
806 if(phi<0.0)
807 phi+=fourier_interval;
808
809 // Add to roots
810 muval.push_back(phi);
811 }
812
813 // Sort the roots
814 std::sort(muval.begin(),muval.end());
815
816 // Sanity check
817 if(!muval.size()) {
818 // Failed. Use polynomial step instead.
819 printf("No root found, falling back to polynomial step.\n");
820 polynomial_step_df(f);
821 return;
822 }
823
824 if( fabs(windowed(0))<sqrt(DBL_EPSILON)) {
825 // Failed. Use polynomial step instead.
826 printf("Derivative value %e %e too small, falling back to polynomial step.\n",fp(0),windowed(0));
827 polynomial_step_df(f);
828 return;
829 }
830
831 // Figure out where the function goes to the wanted direction
832 double findJ;
833 findJ=arma::max(f->getsign()*fv);
834
835 // and the corresponding value of mu is
836 double findmu=mu(0);
837 for(int i=0;i<fourier_length;i++)
838 if(f->getsign()*fv(i)==findJ) {
839 findmu=mu(i);
840 // Stop at closest extremum
841 break;
842 }
843
844 // Find closest value of mu
845 size_t rootind=0;
846 double diffmu=fabs(muval[0]-findmu);
847 for(size_t i=1;i<muval.size();i++)
848 if(fabs(muval[i]-findmu)<diffmu) {
849 rootind=i;
850 diffmu=fabs(muval[i]-findmu);
851 }
852
853 // Optimized step size is
854 double step=muval[rootind];
855
856 // Is the step length in the allowed region?
857 if(step>0.0 && step <=fourier_interval) {
858 // Yes. Calculate the new value
859 arma::cx_mat Wtr=get_rotation(step*f->getsign())*W;
860 UnitaryFunction *newf=f->copy();
861
862 double J=f->getf();
863 double Jtr=newf->cost_func(Wtr);
864
865 // Accept the step?
866 if( f->getsign()*(Jtr-J) > 0.0) {
867 // Yes.
868 delete f;
869 f=newf;
870 // Free memory
871 for(int i=0;i<fourier_length;i++)
872 delete fs[i];
873
874 return;
875 } else
876 delete newf;
877 }
878
879 // If we are still here, then just get the minimum value
880 if(step==0.0 || step > Tmu) {
881 double minval=arma::max(f->getsign()*fv);
882 for(size_t i=0;i<mu.n_elem;i++)
883 if(minval==fv(i)) {
884 delete f;
885 f=fs[i];
886 break;
887 }
888 }
889
890 for(int i=0;i<fourier_length;i++)
891 delete fs[i];
892 }
893
bracket(const arma::cx_mat & X,const arma::cx_mat & Y)894 double bracket(const arma::cx_mat & X, const arma::cx_mat & Y) {
895 return 0.5*std::real(arma::trace(arma::trans(X)*Y));
896 }
897
companion_matrix(const arma::cx_vec & c)898 arma::cx_mat companion_matrix(const arma::cx_vec & c) {
899 if(c.size()<=1) {
900 // Dummy return
901 arma::cx_mat dum;
902 return dum;
903 }
904
905 // Form companion matrix
906 size_t N=c.size()-1;
907 if(c(N)==0.0) {
908 ERROR_INFO();
909 throw std::runtime_error("Coefficient of highest term vanishes!\n");
910 }
911
912 arma::cx_mat companion(N,N);
913 companion.zeros();
914
915 // First row - coefficients normalized to that of highest term.
916 for(size_t j=0;j<N;j++)
917 companion(0,j)=-c(N-(j+1))/c(N);
918 // Fill out the unit matrix part
919 for(size_t j=1;j<N;j++)
920 companion(j,j-1)=1.0;
921
922 return companion;
923 }
924
solve_roots_cplx(const arma::vec & a)925 arma::cx_vec solve_roots_cplx(const arma::vec & a) {
926 return solve_roots_cplx(a*COMPLEX1);
927 }
928
solve_roots_cplx(const arma::cx_vec & a)929 arma::cx_vec solve_roots_cplx(const arma::cx_vec & a) {
930 // Find roots of a_0 + a_1*mu + ... + a_(p-1)*mu^(p-1) = 0.
931
932 // Coefficient of highest order term must be nonzero.
933 size_t r=a.size();
934 while(a(r-1)==0.0 && r>=1)
935 r--;
936
937 if(r==1) {
938 // Zeroth degree - no zeros!
939 arma::cx_vec dummy;
940 dummy.zeros(0);
941 return dummy;
942 }
943
944 // Form companion matrix
945 arma::cx_mat comp=companion_matrix(a.subvec(0,r-1));
946
947 // and diagonalize it
948 arma::cx_vec eigval;
949 arma::cx_mat eigvec;
950 arma::eig_gen(eigval, eigvec, comp);
951
952 // Return the sorted roots
953 return arma::sort(eigval);
954 }
955
solve_roots(const arma::vec & a)956 arma::vec solve_roots(const arma::vec & a) {
957 // Solve the roots
958 arma::cx_vec croots=solve_roots_cplx(a);
959
960 // Collect real roots
961 size_t nreal=0;
962 for(size_t i=0;i<croots.n_elem;i++)
963 if(fabs(std::imag(croots[i]))<10*DBL_EPSILON)
964 nreal++;
965
966 // Real roots
967 arma::vec roots(nreal);
968 size_t ir=0;
969 for(size_t i=0;i<croots.n_elem;i++)
970 if(fabs(std::imag(croots(i)))<10*DBL_EPSILON)
971 roots(ir++)=std::real(croots(i));
972
973 // Sort roots
974 roots=arma::sort(roots);
975
976 return roots;
977 }
978
smallest_positive(const arma::vec & a)979 double smallest_positive(const arma::vec & a) {
980 double step=0.0;
981 for(size_t i=0;i<a.size();i++) {
982 // Omit extremely small steps because they might get you stuck.
983 if(a(i)>sqrt(DBL_EPSILON)) {
984 step=a(i);
985 break;
986 }
987 }
988
989 return step;
990 }
991
derivative_coefficients(const arma::vec & c)992 arma::vec derivative_coefficients(const arma::vec & c) {
993 // Coefficients for polynomial expansion of y'
994 arma::vec cder(c.n_elem-1);
995 for(size_t i=1;i<c.n_elem;i++)
996 cder(i-1)=i*c(i);
997
998 return cder;
999 }
1000
fit_polynomial(const arma::vec & x,const arma::vec & y,int deg)1001 arma::vec fit_polynomial(const arma::vec & x, const arma::vec & y, int deg) {
1002 // Fit function to polynomial of order p: y(x) = a_0 + a_1*x + ... + a_(p-1)*x^(p-1)
1003
1004 if(x.n_elem!=y.n_elem) {
1005 ERROR_INFO();
1006 throw std::runtime_error("x and y have different dimensions!\n");
1007 }
1008 size_t N=x.n_elem;
1009
1010 // Check degree
1011 if(deg<0)
1012 deg=(int) x.size();
1013 if(deg>(int) N) {
1014 ERROR_INFO();
1015 throw std::runtime_error("Underdetermined polynomial!\n");
1016 }
1017
1018 // Form mu matrix
1019 arma::mat mumat(N,deg);
1020 mumat.zeros();
1021 for(size_t i=0;i<N;i++)
1022 for(int j=0;j<deg;j++)
1023 mumat(i,j)=pow(x(i),j);
1024
1025 // Solve for coefficients: mumat * c = y
1026 arma::vec c;
1027 bool solveok=arma::solve(c,mumat,y);
1028 if(!solveok) {
1029 arma::trans(x).print("x");
1030 arma::trans(y).print("y");
1031 mumat.print("Mu");
1032 throw std::runtime_error("Error solving for coefficients a.\n");
1033 }
1034
1035 return c;
1036 }
1037
fit_polynomial_fdf(const arma::vec & x,const arma::vec & y,const arma::vec & dy,int deg)1038 arma::vec fit_polynomial_fdf(const arma::vec & x, const arma::vec & y, const arma::vec & dy, int deg) {
1039 // Fit function and its derivative to polynomial of order p:
1040 // y(x) = a_0 + a_1*x + ... + a_(p-1)*x^(p-1)
1041 // y'(x) = a_1 + ... + (p-1)*a_(p-1)*x^(p-2)
1042 // return coefficients of y'
1043
1044 if(x.n_elem!=y.n_elem) {
1045 ERROR_INFO();
1046 throw std::runtime_error("x and y have different dimensions!\n");
1047 }
1048 if(y.n_elem!=dy.n_elem) {
1049 ERROR_INFO();
1050 throw std::runtime_error("y and dy have different dimensions!\n");
1051 }
1052
1053 // Length of vectors is
1054 size_t N=x.n_elem;
1055 if(deg<0) {
1056 // Default degree of polynomial is
1057 deg=(int) 2*N;
1058 } else {
1059 // We need one more degree so that the derivative is of order deg
1060 deg++;
1061 }
1062
1063 if(deg>(int) (2*N)) {
1064 ERROR_INFO();
1065 throw std::runtime_error("Underdetermined polynomial!\n");
1066 }
1067
1068 // Form mu matrix.
1069 arma::mat mumat(2*N,deg);
1070 mumat.zeros();
1071 // First y(x)
1072 for(size_t i=0;i<N;i++)
1073 for(int j=0;j<deg;j++)
1074 mumat(i,j)=pow(x(i),j);
1075 // Then y'(x)
1076 for(size_t i=0;i<N;i++)
1077 for(int j=1;j<deg;j++)
1078 mumat(i+N,j)=j*pow(x(i),j-1);
1079
1080 // Form rhs vector
1081 arma::vec data(2*N);
1082 data.subvec(0,N-1)=y;
1083 data.subvec(N,2*N-1)=dy;
1084
1085 // Solve for coefficients: mumat * c = data
1086 arma::vec c;
1087 bool solveok=arma::solve(c,mumat,data);
1088 if(!solveok) {
1089 arma::trans(x).print("x");
1090 arma::trans(y).print("y");
1091 arma::trans(dy).print("dy");
1092 mumat.print("Mu");
1093 throw std::runtime_error("Error solving for coefficients a.\n");
1094 }
1095
1096 return c;
1097 }
1098
1099
Brockett(size_t N,unsigned long int seed)1100 Brockett::Brockett(size_t N, unsigned long int seed) : UnitaryFunction(2, true) {
1101 // Get random complex matrix
1102 sigma=randn_mat(N,N,seed)+COMPLEXI*randn_mat(N,N,seed+1);
1103 // Hermitize it
1104 sigma=sigma+arma::trans(sigma);
1105 // Get N matrix
1106 Nmat.zeros(N,N);
1107 for(size_t i=0;i<N;i++)
1108 Nmat(i,i)=i+1;
1109 }
1110
~Brockett()1111 Brockett::~Brockett() {
1112 }
1113
copy() const1114 Brockett* Brockett::copy() const {
1115 return new Brockett(*this);
1116 }
1117
cost_func(const arma::cx_mat & Wv)1118 double Brockett::cost_func(const arma::cx_mat & Wv) {
1119 W=Wv;
1120 f=std::real(arma::trace(arma::trans(W)*sigma*W*Nmat));
1121 return f;
1122 }
1123
cost_der(const arma::cx_mat & Wv)1124 arma::cx_mat Brockett::cost_der(const arma::cx_mat & Wv) {
1125 W=Wv;
1126 return sigma*W*Nmat;
1127 }
1128
cost_func_der(const arma::cx_mat & Wv,double & fv,arma::cx_mat & der)1129 void Brockett::cost_func_der(const arma::cx_mat & Wv, double & fv, arma::cx_mat & der) {
1130 fv=cost_func(Wv);
1131 der=cost_der(Wv);
1132 }
1133
legend() const1134 std::string Brockett::legend() const {
1135 char stat[1024];
1136 sprintf(stat,"%13s %13s", "diag", "unit");
1137 return std::string(stat);
1138 }
1139
status(bool lfmt)1140 std::string Brockett::status(bool lfmt) {
1141 char stat[1024];
1142 if(lfmt)
1143 sprintf(stat,"% .16e % .16e", diagonality(), unitarity());
1144 else
1145 sprintf(stat,"% e % e", diagonality(), unitarity());
1146 return std::string(stat);
1147 }
1148
diagonality() const1149 double Brockett::diagonality() const {
1150 arma::cx_mat WSW=arma::trans(W)*sigma*W;
1151
1152 double off=0.0;
1153 double dg=0.0;
1154
1155 for(size_t i=0;i<WSW.n_cols;i++)
1156 dg+=std::norm(WSW(i,i));
1157
1158 for(size_t i=0;i<WSW.n_cols;i++) {
1159 for(size_t j=0;j<i;j++)
1160 off+=std::norm(WSW(i,j));
1161 for(size_t j=i+1;j<WSW.n_cols;j++)
1162 off+=std::norm(WSW(i,j));
1163 }
1164
1165 return 10*log10(off/dg);
1166 }
1167
unitarity() const1168 double Brockett::unitarity() const {
1169 arma::cx_mat U=W*arma::trans(W);
1170 arma::cx_mat eye(W);
1171 eye.eye();
1172
1173 double norm=pow(arma::norm(U-eye,"fro"),2);
1174 return 10*log10(norm);
1175 }
1176