/** @(#) doglegstep.cpp */
// #include "ap.h"
/** @(#) .h */
/**
* @author Daniel Alba Cuellar
* @author Omar Posada Villarreal
* @version 1.0, /03/2002
*/
#ifndef DOGLEGSTEP_H
#define DOGLEGSTEP_H
//----------------------------------------------------------------------------
#include "../matrixdouble/squaredmatrix.h"
#include "../matrixdouble/vectordouble.h"
#include "../mathpos/mathpos.h"
#include <iomanip>
#include <iostream>
#include <string>
#include <stdexcept> // runtime_exception
using namespace std;
/** Busca una solucion aproximada a
* min g^T s + 1/2 s^T L L^T s
* s E R^n
* sujeto a ||Dx s||2 <= delta
* mediante la seleccion de s en la curca del double dogleg
* descrito abajo de tal forma que ||Dx s||2 = delta,
* o s = SN si Newtlen <= delta. (Dx = diag ((Sx)1, ..., (sx)n).
* Logica uno.
* @author Omar Posada Villarreal, Daniel Alba Cuellar
* @param n Taman~o de vectores y matrices. No se checan dimensiones.
* @param g Gradiente
* @param L Triangular inferior. Se usa para construir el inverso de la Hessiana.
* @param sN = - (L L^T)^-1 g)
* @param Sx
* @param Newtlen = ||Dx sN||2
* @param maxstep Maximo taman~o de paso. Dado por QuasiNewton
* @param delta in/out
* @param firstdog in/out
* false: Cauchylen, eta, sSD, v
* tendran valores actuales en in/out.
* true: Los parametros mencionados no se tomaran en cuenta.
* Inicializa los parametros.
* @param Cauchylen in/out
* @param eta in/out gamma <= eta <=1
* @param sSD in/out
* @param v in/out
* @param s return Paso escogido
* @param Newttaken return
* true: Se usa el paso de Newton
* false: se usa el dogleg step
* @version 1.0, 20/04/2002, sig vers LowerTriangularMatrix &L
* @see Dennis. A6.4.4. DOGSTEP pagina 336
*/
// Ver codigo de preuba al final
void doubleDogLegStep(int n, VectorDouble &g, SquaredMatrix &L,
VectorDouble &sN, VectorDouble &Sx, double Newtlen,
double maxstep,
double &delta, bool &firstdog, double &Cauchylen,
double &eta, VectorDouble &sSD, VectorDouble &v,
VectorDouble &s, bool &Newttaken,
int logType = SHORT_LOG, int waitType = NO_WAIT)
throw (invalid_argument) {
LLOG << "\n{ doubleDogLegStep";
// Es una matriz cuadrada, pero se almacenara como un vector a Sx
VectorDouble Dx(Sx);
VectorDouble DxInv(n, 1);
VectorDouble DxInv_g(n, 1);
VectorDouble Dx_sN(n, 1);
VectorDouble DxInv_sSD(n, 1);
VectorDouble eta_Dx_sN(n, 1);
VectorDouble lambda_v(n, 1);
VectorDouble sSD_lambda_v(n, 1);
double temp, tempv;
double gT_sN;
double alfa, beta, lambda, gamma;
int i, j;
// A6.4.4
if (Newtlen <= delta) { // 1.
// 1T. s es el paso de Newton
Newttaken = true; // 1T.1.
s = sN; // 1T.2. s = sN. antes copyFrom()
delta = Newtlen; // 1T.3.
// A6.4.4. return desde aqui
} else { // 1E.
// el paso de Newton es muy largo, s en curva double dogleg
Newttaken = false; // 1E.1.
if (firstdog) { // 1E.2.
// Calcula la curva double dogleg
// Calcula los 4 parametros de in/out:
// Cauchylen, eta, sSD, v
firstdog = false; // 1E.2.1.
// alfa = ||Dx^-1g||2 ^2
// Dx = diag((Sx)1, ... ,(Sx)n)
// se almacena como Sx E R^n..
// Para implementacion sugerida de expresiones
// que conciernen a Dx, ver seccion 1.5 del apendice
DxInv = Dx.inverse(); // 1E.2.2.
DxInv_g = DxInv.diagonalProduct(g);
//LLOG << "\nDx = " << Dx;
//LLOG << "\nDxInv = " << DxInv;
//LLOG << "\nDxInv_g = " << DxInv_g;
alfa = square( DxInv_g.norm2() );
// 1E.2.3-4. Implementan beta = ||L^T Dx^-2 g||2 ^2
beta = 0.0; // 1E.2.3.
for (i = 1; i <= n; ++i) { // 1E.2.4.
// 1E.2.4.1.
temp = 0.0; // = sum
for (j = 1; j <= n; ++j) {
// Checar vs [][]
temp += (L(j,i) * g(j))
/ square( Sx(j) );
//LLOG << "\nL(j,i) = " << L(j,i);
//LLOG << "\ng(j) = " << g(j);
}
beta += square(temp); // 1E.2.4.2.
}
// ^sSD es el paso de Cauchy en metrica escalada
// ssD en descripcion es Dx^-1 ^sSD
sSD = DxInv_g * (-alfa/beta); // 1E.2.5.
// Cauchylen = ||^sSD||2
Cauchylen = pow(alfa, 1.5) / beta; // 1E.2.6.
LLOG << "\ndelta = " << delta; // didac
LLOG << "\n||s^CP||2. Cauchylen = " << Cauchylen;
// didac
LLOG << "\nDentro de la region de confianza?"
<< " CP < delta = "
<< toStr(Cauchylen < delta);
// didac
// Son vectores g.transpose().multiply(sN);
gT_sN = g * sN; // innerProduct()
// eta <= 1. ver seccion 6.4.2.
// 1E.2.7.
gamma = square(alfa) / (beta * fabs(gT_sN)); // didac
eta = 0.2 + (0.8 * gamma);
LLOG << "\ngamma <= eta <= 1. gamma = " << gamma;
// didac
// 1E.2.8.
// ^v es (eta sN - sSD) en la metrica escalada.
// ^v = eta Dx sN - sSD
Dx_sN = Dx.diagonalProduct(sN); // No usar Dx*sN
eta_Dx_sN = Dx_sN * eta;
v = eta_Dx_sN - sSD;
if ( areAlmostEqual(delta, -1) ) { // 1E.2.9.
// 1a iteracion, y el usuario no dio una region
// confiable inicial
delta = min(Cauchylen, maxstep);
}
} // 1E.2.
if ((eta * Newtlen) <= delta) { // 1.E.3a.
// Toma un paso intermedio en la direccion de Newton
s = sN * (delta / Newtlen); // 1.E.3a.1.
// A6.4.4 regresa de aqui
} else if (Cauchylen >= delta) { // 1.E.3b
// Toma un paso en la direccion descendiente
// de mayor profundidad
// 1.E.3b.1
// s = (delta / Cauchylen) * Dx^-1 ^sSD
DxInv = Dx.inverse();
DxInv_sSD = Dx.diagonalProduct(sSD);
s = DxInv_sSD * (delta / Cauchylen);
// A6.4.4 regresa de aqui
} else { // 1.E.3c
// Calcula la combinacion convexa de sSD y eta sN
// (= DX^-1 ^sSD + lambda Dx^-1 ^v)
// que tiene largo escalado = delta
temp = v * sSD; // 1.E.3c.1 innerProduct()
tempv = v * v; // 1.E.3c.2 innerProduct()
// 1.E.3c.3
lambda = (sqrt(
square(temp)
- (tempv * (
square(Cauchylen)
- square(delta))))
- temp) / tempv;
// 1.E.3c.4.
// s = Dx^-1 (^sSD + lambda ^v)
DxInv = Dx.inverse();
lambda_v = v * lambda;
sSD_lambda_v = sSD + lambda_v;
s = DxInv.diagonalProduct(sSD_lambda_v);
} // 1E.3c.
} // 1.E.
LLOG << "\n } doubleDogLegStep";
}
//----------------------------------------------------------------------------
#endif
/* CODIGO DEL MAIN PARA PROBAR TESTDOGLEG
int main() {
try {
int n = 2; //dimension
double fc; // escalar
VectorDouble xc(n, 1);
VectorDouble Sx(n, 1);
SquaredMatrix H(n, 1);
Function1 func;
// Proceso
VectorDouble g(n, 1);
SquaredMatrix L(n, 1);
VectorDouble sN(n, 1);
//VectorDouble Sx(n, 1);
double Newtlen;
double maxstep;
double delta;
bool firstdog;
double Cauchylen;
double eta;
VectorDouble sSD(n, 1);
VectorDouble v(n, 1);
VectorDouble s(n, 1);
bool Newttaken;
SquaredMatrix LT_L(n, 1);
SquaredMatrix LT_LInv(n, 1);
/*
xc = [ 1 1 ]
H(xc) =
[ 14 0 ]
[ 0 2 ]
H^-1(xc) =
[ 0.0714286 0 ]
[ 0 0.5 ]
*/
/* xc[1] = 1;
xc[2] = 1;
Sx[1] = 1;
Sx[2] = 1;
H[1][1] = 14;
H[1][2] = 0;
H[2][1] = 0;
H[2][2] = 2;
if (!factorizeCholesky(n, H, L) ) {
throw runtime_error("No se pudo factorizar.");
}
g[1] = 6.0;
g[2] = 2.0;
//LT_L = L.transpose() * L;
LT_LInv[1][1] = 1.0/14.0;
LT_LInv[2][2] = 1.0/2.0;
//error sN = (LT_LInv * g) * (-1.0);
sN[1] = -0.428571;
sN[2] = -1.0;
// Es una matriz cuadrada, pero se almacenara como un vector a Sx
VectorDouble Dx(Sx);
VectorDouble Dx_sN(n, 1);
Dx_sN = Dx.diagonalProduct(sN);
Newtlen = Dx_sN.norm2();
// Dado por el usuario
maxstep = 1.0;
// Dado por el usuario
delta = 0.75; //ejemplo
firstdog = true; // init updTrusReg
cout << "\nPunto inicial. xc: " << xc;
cout << "\nHessiana. H:\n" << H;
waitEnter();
cout << "\n-----IN-----";
cout << "\nDimension. n = " << n;
cout << "\ngrad(f(xc)). g = " << g;
cout << "\nTriangular inferior del Hessiano. L:\n" << L;
cout << "\nsc^N = -Hc^-1 grad(f(xc)). sN:\n" << sN;
cout << "\nEscala metrica. Sx: " << Sx;
cout << "\nPaso de Newton. Newtlen = " << Newtlen;
cout << "\nMaxima delta. maxstep = " << maxstep;
cout << "\nRadio de la region de confianza. delta = " << delta;
cout << "\ntrue: 1a. iteracion. firstdog = " << firstdog;
cout << "\n-----doubleDogLegStep-----";
doubleDogLegStep(n, g, L,
sN, Sx, Newtlen,
maxstep, delta, firstdog,
Cauchylen, eta, sSD,
v,
s, Newttaken);
cout << "\n-----IN/OUT-----";
cout << "\n||s^CP||2. Cauchylen = " << Cauchylen;
cout << "\nFactor para pasar de CP a ^N. eta = " << eta;
cout << "\ns^CP, vector de xc a CP. sSD = " << sSD;
cout << "\neta s^(^N) - s^CP. Vector para pasar de CP a ^N. v = " << v;
cout << "\n-----OUT-----";
cout << "\ns^CP + lambda v. Paso. s = " << s;
cout << "\ntrue: se tomo el paso de Newton. Newttaken = "
<< toString(Newttaken);
/* Resultado de prueba. xc = (2, 2)T. H(xc):
* [13 -4]
* [-4 1 ]
* Como esta especificado el algoritmo, solo se encuentra la H
* para la diagonal y la matriz triangular superior.
* El resto no se toca porque es simetrica.
*/
/* } catch (exception e) {
cout << "\nmain: Excepcion " << e.what();
}
// Terminacion
cout << "\n***Fin***";
waitEnter();
//si no hay pausa cin.get();
//cin.get();
return 0;
} Fin TEST-----------------------------------------------------------------
*/
// Fin------------------------------------------------------------------------