/** * Dense trust-region unconstrained minimization using Dogleg steps and BFGS to estimate the Hessian. * @param config Trust region configuration * @return The new optimization routine */ public static UnconstrainedMinimization doglegBFGS( @Nullable ConfigTrustRegion config ) { if( config == null ) config = new ConfigTrustRegion(); HessianBFGS hessian = new HessianBFGS_DDRM(true); TrustRegionUpdateDogleg_F64<DMatrixRMaj> update = new TrustRegionUpdateDogleg_F64<>(); UnconMinTrustRegionBFGS_F64 alg = new UnconMinTrustRegionBFGS_F64(update,hessian); alg.configure(config); return alg; }
protected void combinedStep(double regionRadius, DMatrixRMaj step) { // find the Cauchy point CommonOps_DDRM.scale(-distanceCauchy, direction, stepCauchy); stepLength = regionRadius; // touches the trust region double distancePtoGN = SpecializedOps_DDRM.diffNormF(stepCauchy,stepGN); double f = fractionCauchyToGN(distanceCauchy,distanceGN,distancePtoGN,regionRadius); CommonOps_DDRM.add(1-f,stepCauchy,f,stepGN,step); predictedReduction = owner.computePredictedReduction(step); }
@Override public void computeUpdate(DMatrixRMaj step, double regionRadius) { if( positiveDefinite ) { // If the GN solution is inside the trust region it should use that solution if( distanceGN <= regionRadius ) { if( verbose != null ) verbose.println(" newton"); gaussNewtonStep(step); } else if( distanceCauchy >= regionRadius ) { if( verbose != null ) verbose.println(" cauchy"); // if the trust region comes before the Cauchy point then perform the cauchy step cauchyStep(regionRadius, step); } else { if( verbose != null ) verbose.println(" combined"); // the solution lies on the line connecting Cauchy and GN combinedStep(regionRadius, step); } } else { if( verbose != null ) verbose.println(" not positive-definite. gBg="+gBg); // Cauchy step for negative semi-definite systems stepLength = regionRadius; CommonOps_DDRM.scale(-stepLength, direction,step); predictedReduction = stepLength*owner.gradientNorm - 0.5*stepLength*stepLength*gBg; } }
@Test public void computeUpdate_cauchy_after() { TrustRegionUpdateDogleg_F64<DMatrixRMaj> alg = new TrustRegionUpdateDogleg_F64<>(); MockTrustRegionBase owner = new MockTrustRegionBase(alg); alg.initialize(owner,2,0); owner.hessian().set(new double[][]{{2,0.1},{0.1,1.5}}); // have Gn and cauchy lie along a line so the math is easy owner.gradientNorm = 2.1; setGradient(alg.owner.gradient,-1,0,owner.gradientNorm); CommonOps_DDRM.divide(owner.gradient,owner.gradientNorm,alg.direction); alg.gBg = owner.hessian.innerVectorHessian(alg.direction); alg.positiveDefinite = true; alg.distanceGN = 5; alg.distanceCauchy = 3; alg.stepCauchy.set(new double[][]{{3},{0}}); alg.stepGN.set(new double[][]{{5},{0}}); double radius = 2; DMatrixRMaj p = new DMatrixRMaj(2,1); alg.computeUpdate(p,radius); assertEquals(owner.computePredictedReduction(p),alg.getPredictedReduction(),UtilEjml.TEST_F64); assertEquals(radius, alg.getStepLength(), UtilEjml.TEST_F64); assertEquals(radius, NormOps_DDRM.normF(p), UtilEjml.TEST_F64); }
/** * See if it correctly identifies a non SPD matrix */ @Test public void initializeUpdate_spd() { TrustRegionUpdateDogleg_F64<DMatrixRMaj> alg = new TrustRegionUpdateDogleg_F64<>(); MockTrustRegionBase owner = new MockTrustRegionBase(alg); alg.initialize(owner,2,0); owner.hessian().set(new double[][]{{1,0},{0,1}}); owner.gradientNorm = 1; owner.gradient.set(new double[][]{{1},{0}}); alg.initializeUpdate(); assertTrue(alg.positiveDefinite); owner.hessian().set(new double[][]{{0,1},{1,0}}); alg.initializeUpdate(); assertFalse(alg.positiveDefinite); }
@Override public void initializeUpdate() { // Scale the gradient vector to make it less likely to overflow/underflow CommonOps_DDRM.divide(owner.gradient,owner.gradientNorm, direction); gBg = owner.hessian.innerVectorHessian(direction); if(UtilEjml.isUncountable(gBg)) throw new OptimizationException("Uncountable. gBg="+gBg); // see if it's positive definite and Gauss-Newton can be solved if( gBg > 0 && solveGaussNewtonPoint(stepGN) ) { positiveDefinite = true; // length of the Cauchy step when computed without constraints distanceCauchy = owner.gradientNorm/gBg; // p_gn = -inv(B)*g CommonOps_DDRM.scale(-1, stepGN); distanceGN = NormOps_DDRM.normF(stepGN); } else { positiveDefinite = false; } }
@Test public void computeUpdate_NegativeDefinite() { TrustRegionUpdateDogleg_F64<DMatrixRMaj> alg = new TrustRegionUpdateDogleg_F64<>(); MockTrustRegionBase owner = new MockTrustRegionBase(alg); alg.initialize(owner,2,0); double radius = 2; setGradient(alg.owner.gradient,-1,0,1000); alg.owner.gradientNorm = NormOps_DDRM.normF(alg.owner.gradient); CommonOps_DDRM.divide(owner.gradient,owner.gradientNorm,alg.direction); owner.hessian().set(new double[][]{{-2,0.1},{0.1,-1.5}}); alg.gBg = owner.hessian.innerVectorHessian(alg.direction); alg.positiveDefinite = false; alg.owner.fx = 1000; DMatrixRMaj p = new DMatrixRMaj(2,1); alg.computeUpdate(p,radius); assertEquals(2,p.get(0,0), UtilEjml.TEST_F64); assertEquals(0,p.get(1,0), UtilEjml.TEST_F64); assertEquals(owner.computePredictedReduction(p),alg.getPredictedReduction(),UtilEjml.TEST_F64); assertEquals(radius, alg.getStepLength(), UtilEjml.TEST_F64); assertEquals(radius, NormOps_DDRM.normF(p), UtilEjml.TEST_F64); }
@Override public void initializeUpdate() { // Scale the gradient vector to make it less likely to overflow/underflow CommonOps_DDRM.divide(owner.gradient,owner.gradientNorm, direction); gBg = owner.hessian.innerVectorHessian(direction); if(UtilEjml.isUncountable(gBg)) throw new OptimizationException("Uncountable. gBg="+gBg); // see if it's positive definite and Gauss-Newton can be solved if( gBg > 0 && solveGaussNewtonPoint(stepGN) ) { positiveDefinite = true; // length of the Cauchy step when computed without constraints distanceCauchy = owner.gradientNorm/gBg; // p_gn = -inv(B)*g CommonOps_DDRM.scale(-1, stepGN); distanceGN = NormOps_DDRM.normF(stepGN); } else { positiveDefinite = false; } }
/** * The Gauss-Newton step is inside the trust region so we will just use that */ @Test public void computeUpdate_gn_inside() { TrustRegionUpdateDogleg_F64<DMatrixRMaj> alg = new TrustRegionUpdateDogleg_F64<>(); MockTrustRegionBase owner = new MockTrustRegionBase(alg); alg.initialize(owner,2,0); alg.positiveDefinite = true; owner.hessian().set(new double[][]{{2,0.1},{0.1,1.5}}); setGradient(alg.owner.gradient,-1,0,1.5); alg.owner.gradientNorm = NormOps_DDRM.normF(alg.owner.gradient); alg.gBg = owner.hessian.innerVectorHessian(alg.direction); alg.stepGN.set(new double[][]{{1},{2}}); alg.distanceGN = NormOps_DDRM.normF(alg.stepGN); DMatrixRMaj p = new DMatrixRMaj(2,1); alg.computeUpdate(p,4); assertTrue(MatrixFeatures_DDRM.isIdentical(p,alg.stepGN,UtilEjml.TEST_F64)); assertEquals(owner.computePredictedReduction(p),alg.getPredictedReduction(),UtilEjml.TEST_F64); assertEquals(alg.distanceGN, alg.getStepLength(), UtilEjml.TEST_F64); assertEquals(alg.distanceGN, NormOps_DDRM.normF(p), UtilEjml.TEST_F64); }
/** * Dense trust-region unconstrained minimization using Dogleg steps and BFGS to estimate the Hessian. * @param config Trust region configuration * @return The new optimization routine */ public static UnconstrainedMinimization doglegBFGS( @Nullable ConfigTrustRegion config ) { if( config == null ) config = new ConfigTrustRegion(); HessianBFGS hessian = new HessianBFGS_DDRM(true); TrustRegionUpdateDogleg_F64<DMatrixRMaj> update = new TrustRegionUpdateDogleg_F64<>(); UnconMinTrustRegionBFGS_F64 alg = new UnconMinTrustRegionBFGS_F64(update,hessian); alg.configure(config); return alg; }
@Override public void computeUpdate(DMatrixRMaj step, double regionRadius) { if( positiveDefinite ) { // If the GN solution is inside the trust region it should use that solution if( distanceGN <= regionRadius ) { if( verbose != null ) verbose.println(" newton"); gaussNewtonStep(step); } else if( distanceCauchy >= regionRadius ) { if( verbose != null ) verbose.println(" cauchy"); // if the trust region comes before the Cauchy point then perform the cauchy step cauchyStep(regionRadius, step); } else { if( verbose != null ) verbose.println(" combined"); // the solution lies on the line connecting Cauchy and GN combinedStep(regionRadius, step); } } else { if( verbose != null ) verbose.println(" not positive-definite. gBg="+gBg); // Cauchy step for negative semi-definite systems stepLength = regionRadius; CommonOps_DDRM.scale(-stepLength, direction,step); predictedReduction = stepLength*owner.gradientNorm - 0.5*stepLength*stepLength*gBg; } }
protected void combinedStep(double regionRadius, DMatrixRMaj step) { // find the Cauchy point CommonOps_DDRM.scale(-distanceCauchy, direction, stepCauchy); stepLength = regionRadius; // touches the trust region double distancePtoGN = SpecializedOps_DDRM.diffNormF(stepCauchy,stepGN); double f = fractionCauchyToGN(distanceCauchy,distanceGN,distancePtoGN,regionRadius); CommonOps_DDRM.add(1-f,stepCauchy,f,stepGN,step); predictedReduction = owner.computePredictedReduction(step); }
@Test public void computeUpdate_cauchy_inside() { TrustRegionUpdateDogleg_F64<DMatrixRMaj> alg = new TrustRegionUpdateDogleg_F64<>(); MockTrustRegionBase owner = new MockTrustRegionBase(alg); alg.initialize(owner,2,0); // have Gn and cauchy lie along a line so the math is easy double radius = 2; owner.hessian().set(new double[][]{{2,0.1},{0.1,1.5}}); alg.owner.gradientNorm = 1.5; setGradient(alg.owner.gradient,-1,0,alg.owner.gradientNorm); CommonOps_DDRM.divide(owner.gradient,owner.gradientNorm,alg.direction); alg.gBg = owner.hessian.innerVectorHessian(alg.direction); alg.positiveDefinite = true; alg.stepGN.set(new double[][]{{5},{0}}); alg.distanceGN = 5; alg.stepCauchy.set(new double[][]{{1.5},{0}}); alg.distanceCauchy = NormOps_DDRM.normF(alg.stepCauchy); alg.gBg = 1; DMatrixRMaj p = new DMatrixRMaj(2,1); alg.computeUpdate(p,radius); assertEquals(2,p.get(0,0), UtilEjml.TEST_F64); assertEquals(0,p.get(1,0), UtilEjml.TEST_F64); assertEquals(owner.computePredictedReduction(p),alg.getPredictedReduction(),UtilEjml.TEST_F64); assertEquals(radius, alg.getStepLength(), UtilEjml.TEST_F64); assertEquals(radius, NormOps_DDRM.normF(p), UtilEjml.TEST_F64); }
/** * Creates a sparse Schur Complement trust region optimization using dogleg steps. * * @see UnconLeastSqTrustRegionSchur_F64 * * @param config Trust region configuration * @return The new optimization routine */ public static UnconstrainedLeastSquaresSchur<DMatrixSparseCSC> doglegSchur( @Nullable ConfigTrustRegion config ) { if( config == null ) config = new ConfigTrustRegion(); HessianSchurComplement_DSCC hessian = new HessianSchurComplement_DSCC(); TrustRegionUpdateDogleg_F64<DMatrixSparseCSC> update = new TrustRegionUpdateDogleg_F64<>(); UnconLeastSqTrustRegionSchur_F64<DMatrixSparseCSC> alg = new UnconLeastSqTrustRegionSchur_F64<>(update,hessian); alg.configure(config); return alg; }
/** * Easy to derive solutions */ @Test public void fractionCauchyToGN_easy() { // Everything lies along a line double lengthPtoGN = 2; double found = TrustRegionUpdateDogleg_F64.fractionCauchyToGN(2,4,lengthPtoGN,2.5); assertEquals(0.5/lengthPtoGN,found,UtilEjml.TEST_F64); }
/** * Creates a sparse Schur Complement trust region optimization using dogleg steps. * * @see UnconLeastSqTrustRegionSchur_F64 * * @param config Trust region configuration * @return The new optimization routine */ public static UnconstrainedLeastSquaresSchur<DMatrixSparseCSC> doglegSchur( @Nullable ConfigTrustRegion config ) { if( config == null ) config = new ConfigTrustRegion(); HessianSchurComplement_DSCC hessian = new HessianSchurComplement_DSCC(); TrustRegionUpdateDogleg_F64<DMatrixSparseCSC> update = new TrustRegionUpdateDogleg_F64<>(); UnconLeastSqTrustRegionSchur_F64<DMatrixSparseCSC> alg = new UnconLeastSqTrustRegionSchur_F64<>(update,hessian); alg.configure(config); return alg; }
/** * Randomly generate points in 2D and circles. Then see if a valid length can be found */ @Test public void fractionCauchyToGN_random() { for (int i = 0; i < 200; i++) { double r = rand.nextDouble()+1; double lengthP = 0.01+rand.nextDouble()*0.99*r; double lengthGN = r + rand.nextDouble(); double angleP = rand.nextDouble()*Math.PI*2.0; double angleGN = rand.nextDouble()*Math.PI*2.0; double x_p = Math.cos(angleP)*lengthP; double y_p = Math.sin(angleP)*lengthP; double x_gn = Math.cos(angleGN)*lengthGN; double y_gn = Math.sin(angleGN)*lengthGN; double dx = x_gn-x_p; double dy = y_gn-y_p; double lengthPtoGN = Math.sqrt(dx*dx + dy*dy); double fraction = TrustRegionUpdateDogleg_F64.fractionCauchyToGN(lengthP,lengthGN,lengthPtoGN,r); double x = x_p + fraction*dx; double y = y_p + fraction*dy; double found = Math.sqrt(x*x + y*y); assertEquals(r,found, UtilEjml.TEST_F64); } }
protected UnconMinTrustRegionBFGS_F64 createAlg() { TrustRegionUpdateDogleg_F64 dogleg = new TrustRegionUpdateDogleg_F64(); HessianBFGS_DDRM hessian = new HessianBFGS_DDRM(true); return new UnconMinTrustRegionBFGS_F64(dogleg,hessian); } }
/** * Creates a sparse trust region optimization using dogleg steps. * * @see UnconLeastSqTrustRegion_F64 * * @param config Trust region configuration * @return The new optimization routine */ public static UnconstrainedLeastSquares<DMatrixSparseCSC> dogleg( @Nullable ConfigTrustRegion config) { if( config == null ) config = new ConfigTrustRegion(); LinearSolverSparse<DMatrixSparseCSC,DMatrixRMaj> solver = LinearSolverFactory_DSCC.cholesky(FillReducing.NONE); HessianLeastSquares_DSCC hessian = new HessianLeastSquares_DSCC(solver); MatrixMath_DSCC math = new MatrixMath_DSCC(); TrustRegionUpdateDogleg_F64<DMatrixSparseCSC> update = new TrustRegionUpdateDogleg_F64<>(); UnconLeastSqTrustRegion_F64<DMatrixSparseCSC> alg = new UnconLeastSqTrustRegion_F64<>(update,hessian,math); alg.configure(config); return alg; }
/** * Creates a sparse trust region optimization using dogleg steps. * * @see UnconLeastSqTrustRegion_F64 * * @param config Trust region configuration * @return The new optimization routine */ public static UnconstrainedLeastSquares<DMatrixSparseCSC> dogleg( @Nullable ConfigTrustRegion config) { if( config == null ) config = new ConfigTrustRegion(); LinearSolverSparse<DMatrixSparseCSC,DMatrixRMaj> solver = LinearSolverFactory_DSCC.cholesky(FillReducing.NONE); HessianLeastSquares_DSCC hessian = new HessianLeastSquares_DSCC(solver); MatrixMath_DSCC math = new MatrixMath_DSCC(); TrustRegionUpdateDogleg_F64<DMatrixSparseCSC> update = new TrustRegionUpdateDogleg_F64<>(); UnconLeastSqTrustRegion_F64<DMatrixSparseCSC> alg = new UnconLeastSqTrustRegion_F64<>(update,hessian,math); alg.configure(config); return alg; }