-
Notifications
You must be signed in to change notification settings - Fork 43
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
Showing
60 changed files
with
5,208 additions
and
45 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
135 changes: 135 additions & 0 deletions
135
...main/java/org/hipparchus/optim/nonlinear/vector/constrained/ADMMQPConvergenceChecker.java
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,135 @@ | ||
| /* | ||
| * Licensed to the Hipparchus project under one or more | ||
| * contributor license agreements. See the NOTICE file distributed with | ||
| * this work for additional information regarding copyright ownership. | ||
| * The Hipparchus project licenses this file to You under the Apache License, Version 2.0 | ||
| * (the "License"); you may not use this file except in compliance with | ||
| * the License. You may obtain a copy of the License at | ||
| * | ||
| * https://www.apache.org/licenses/LICENSE-2.0 | ||
| * | ||
| * Unless required by applicable law or agreed to in writing, software | ||
| * distributed under the License is distributed on an "AS IS" BASIS, | ||
| * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. | ||
| * See the License for the specific language governing permissions and | ||
| * limitations under the License. | ||
| */ | ||
| package org.hipparchus.optim.nonlinear.vector.constrained; | ||
|
|
||
| import org.hipparchus.linear.RealMatrix; | ||
| import org.hipparchus.linear.RealVector; | ||
| import org.hipparchus.optim.ConvergenceChecker; | ||
| import org.hipparchus.optim.OptimizationData; | ||
| import org.hipparchus.util.FastMath; | ||
|
|
||
| /** Convergence Checker for ADMM QP Optimizer. | ||
| * @since 3.1 | ||
| */ | ||
| public class ADMMQPConvergenceChecker implements ConvergenceChecker<LagrangeSolution>, OptimizationData { | ||
|
|
||
| /** Quadratic term matrix. */ | ||
| private final RealMatrix h; | ||
|
|
||
| /** Constraint coefficients matrix. */ | ||
| private final RealMatrix a; | ||
|
|
||
| /** Linear term matrix. */ | ||
| private final RealVector q; | ||
|
|
||
| /** Absolute tolerance for convergence. */ | ||
| private final double epsAbs; | ||
|
|
||
| /** Relative tolerance for convergence. */ | ||
| private final double epsRel; | ||
|
|
||
| /** Convergence indicator. */ | ||
| private boolean converged; | ||
|
|
||
| /** Simple constructor. | ||
| * @param h quadratic term matrix | ||
| * @param a constraint coefficients matrix | ||
| * @param q linear term matrix | ||
| * @param epsAbs | ||
| * @param epsRel | ||
| */ | ||
| ADMMQPConvergenceChecker(final RealMatrix h, final RealMatrix a, final RealVector q, | ||
| final double epsAbs, final double epsRel) { | ||
| this.h = h; | ||
| this.a = a; | ||
| this.q = q; | ||
| this.epsAbs = epsAbs; | ||
| this.epsRel = epsRel; | ||
| this.converged = false; | ||
| } | ||
|
|
||
| /** {@inheritDoc} */ | ||
| @Override | ||
| public boolean converged(final int i, final LagrangeSolution previous, final LagrangeSolution current) { | ||
| return converged; | ||
| } | ||
|
|
||
| /** Evaluate convergence. | ||
| * @param rp primal residual | ||
| * @param rd dual residual | ||
| * @param maxPrimal primal vectors max | ||
| * @param maxDual dual vectors max | ||
| * @return true of convergence has been reached | ||
| */ | ||
| public boolean converged(final double rp, final double rd, final double maxPrimal, final double maxDual) { | ||
| boolean result = false; | ||
|
|
||
| if (rp <= epsPrimalDual(maxPrimal) && rd <= epsPrimalDual(maxDual)) { | ||
| result = true; | ||
| converged = true; | ||
| } | ||
| return result; | ||
| } | ||
|
|
||
| /** Compute primal residual. | ||
| * @param x primal problem solution | ||
| * @param z auxiliary variable | ||
| * @return primal residual | ||
| */ | ||
| public double residualPrime(final RealVector x, final RealVector z) { | ||
| return a.operate(x).subtract(z).getLInfNorm(); | ||
| } | ||
|
|
||
| /** Compute dual residual. | ||
| * @param x primal problem solution | ||
| * @param y dual problem solution | ||
| * @return dual residual | ||
| */ | ||
| public double residualDual(final RealVector x, final RealVector y) { | ||
| return q.add(a.transpose().operate(y)).add(h.operate(x)).getLInfNorm(); | ||
| } | ||
|
|
||
| /** Compute primal vectors max. | ||
| * @param x primal problem solution | ||
| * @param z auxiliary variable | ||
| * @return primal vectors max | ||
| */ | ||
| public double maxPrimal(final RealVector x, final RealVector z) { | ||
| return FastMath.max(a.operate(x).getLInfNorm(), z.getLInfNorm()); | ||
| } | ||
|
|
||
| /** Compute dual vectors max. | ||
| * @param x primal problem solution | ||
| * @param y dual problem solution | ||
| * @return dual vectors max | ||
| */ | ||
| public double maxDual(final RealVector x, final RealVector y) { | ||
| return FastMath.max(FastMath.max(h.operate(x).getLInfNorm(), | ||
| a.transpose().operate(y).getLInfNorm()), | ||
| q.getLInfNorm()); | ||
| } | ||
|
|
||
| /** Combine absolute and relative tolerances. | ||
| * @param maxPrimalDual either {@link #maxPrimal(RealVector, RealVector)} | ||
| * or {@link #maxDual(RealVector, RealVector)} | ||
| * @return global tolerance | ||
| */ | ||
| private double epsPrimalDual(final double maxPrimalDual) { | ||
| return epsAbs + epsRel * maxPrimalDual; | ||
| } | ||
|
|
||
| } |
171 changes: 171 additions & 0 deletions
171
...chus-optim/src/main/java/org/hipparchus/optim/nonlinear/vector/constrained/ADMMQPKKT.java
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,171 @@ | ||
| /* | ||
| * Licensed to the Hipparchus project under one or more | ||
| * contributor license agreements. See the NOTICE file distributed with | ||
| * this work for additional information regarding copyright ownership. | ||
| * The Hipparchus project licenses this file to You under the Apache License, Version 2.0 | ||
| * (the "License"); you may not use this file except in compliance with | ||
| * the License. You may obtain a copy of the License at | ||
| * | ||
| * https://www.apache.org/licenses/LICENSE-2.0 | ||
| * | ||
| * Unless required by applicable law or agreed to in writing, software | ||
| * distributed under the License is distributed on an "AS IS" BASIS, | ||
| * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. | ||
| * See the License for the specific language governing permissions and | ||
| * limitations under the License. | ||
| */ | ||
| package org.hipparchus.optim.nonlinear.vector.constrained; | ||
|
|
||
|
|
||
| import org.hipparchus.linear.ArrayRealVector; | ||
| import org.hipparchus.linear.DecompositionSolver; | ||
| import org.hipparchus.linear.EigenDecompositionSymmetric; | ||
| import org.hipparchus.linear.MatrixUtils; | ||
| import org.hipparchus.linear.RealMatrix; | ||
| import org.hipparchus.linear.RealVector; | ||
| import org.hipparchus.util.FastMath; | ||
|
|
||
| /** Alternative Direction Method of Multipliers Solver. | ||
| * @since 3.1 | ||
| */ | ||
| public class ADMMQPKKT implements KarushKuhnTuckerSolver<ADMMQPSolution> { | ||
|
|
||
| /** Square matrix of weights for quadratic terms. */ | ||
| private RealMatrix H; | ||
|
|
||
| /** Vector of weights for linear terms. */ | ||
| private RealVector q; | ||
|
|
||
| /** constraints coefficients matrix. */ | ||
| private RealMatrix A; | ||
|
|
||
| /** Regularization term sigma for Karush–Kuhn–Tucker solver. */ | ||
| private double sigma; | ||
|
|
||
| /** TBC. */ | ||
| private RealMatrix R; | ||
|
|
||
| /** Inverse of R. */ | ||
| private RealMatrix Rinv; | ||
|
|
||
| /** Lower bound. */ | ||
| private RealVector lb; | ||
|
|
||
| /** Upper bound. */ | ||
| private RealVector ub; | ||
|
|
||
| /** Alpha filter for ADMM iteration. */ | ||
| private double alpha; | ||
|
|
||
| /** Constrained problem KKT matrix. */ | ||
| private RealMatrix M; | ||
|
|
||
| /** Solver for M. */ | ||
| private DecompositionSolver dsX; | ||
|
|
||
| /** Simple constructor. | ||
| * <p> | ||
| * BEWARE, nothing is initialized here, it is {@link #initialize(RealMatrix, RealMatrix, | ||
| * RealVector, int, RealVector, RealVector, double, double, double) initialize} <em>must</em> | ||
| * be called before using the instance. | ||
| * </p> | ||
| */ | ||
| ADMMQPKKT() { | ||
| // nothing initialized yet! | ||
| } | ||
|
|
||
| /** {@inheritDoc} */ | ||
| @Override | ||
| public ADMMQPSolution solve(RealVector b1, final RealVector b2) { | ||
| RealVector z = dsX.solve(new ArrayRealVector((ArrayRealVector) b1,b2)); | ||
| return new ADMMQPSolution(z.getSubVector(0,b1.getDimension()), z.getSubVector(b1.getDimension(), b2.getDimension())); | ||
| } | ||
|
|
||
| public void updateSigmaRho(double newSigma, int me, double rho) { | ||
| this.sigma = newSigma; | ||
| this.H = H.add(MatrixUtils.createRealIdentityMatrix(H.getColumnDimension()).scalarMultiply(newSigma)); | ||
| createPenaltyMatrix(me, rho); | ||
| M = MatrixUtils.createRealMatrix(H.getRowDimension() + A.getRowDimension(), | ||
| H.getRowDimension() + A.getRowDimension()); | ||
| M.setSubMatrix(H.getData(), 0,0); | ||
| M.setSubMatrix(A.getData(), H.getRowDimension(),0); | ||
| M.setSubMatrix(A.transpose().getData(), 0, H.getRowDimension()); | ||
| M.setSubMatrix(Rinv.scalarMultiply(-1.0).getData(), H.getRowDimension(),H.getRowDimension()); | ||
| dsX = new EigenDecompositionSymmetric(M).getSolver(); | ||
| } | ||
|
|
||
| public void initialize(RealMatrix newH, RealMatrix newA, RealVector newQ, | ||
| int me, RealVector newLb, RealVector newUb, | ||
| double rho, double newSigma, double newAlpha) { | ||
| this.lb = newLb; | ||
| this.ub = newUb; | ||
| this.alpha = newAlpha; | ||
| this.sigma = newSigma; | ||
| this.H = newH.add(MatrixUtils.createRealIdentityMatrix(newH.getColumnDimension()).scalarMultiply(newSigma)); | ||
| this.A = newA.copy(); | ||
| this.q = newQ.copy(); | ||
| createPenaltyMatrix(me, rho); | ||
|
|
||
| M = MatrixUtils.createRealMatrix(newH.getRowDimension() + newA.getRowDimension(), | ||
| newH.getRowDimension() + newA.getRowDimension()); | ||
| M.setSubMatrix(newH.getData(),0,0); | ||
| M.setSubMatrix(newA.getData(),newH.getRowDimension(),0); | ||
| M.setSubMatrix(newA.transpose().getData(),0,newH.getRowDimension()); | ||
| M.setSubMatrix(Rinv.scalarMultiply(-1.0).getData(),newH.getRowDimension(),newH.getRowDimension()); | ||
| dsX = new EigenDecompositionSymmetric(M).getSolver(); | ||
| } | ||
|
|
||
| private void createPenaltyMatrix(int me, double rho) { | ||
| this.R = MatrixUtils.createRealIdentityMatrix(A.getRowDimension()); | ||
|
|
||
| for (int i = 0; i < R.getRowDimension(); i++) { | ||
| if (i < me) { | ||
| R.setEntry(i, i, rho * 1000.0); | ||
|
|
||
| } else { | ||
| R.setEntry(i, i, rho); | ||
|
|
||
| } | ||
| } | ||
| this.Rinv = MatrixUtils.inverse(R); | ||
| } | ||
|
|
||
| /** {@inheritDoc} */ | ||
| @Override | ||
| public ADMMQPSolution iterate(RealVector... previousSol) { | ||
| double onealfa = 1.0 - alpha; | ||
| //SAVE OLD VALUE | ||
| RealVector xold = previousSol[0].copy(); | ||
| RealVector yold = previousSol[1].copy(); | ||
| RealVector zold = previousSol[2].copy(); | ||
|
|
||
| //UPDATE RIGHT VECTOR | ||
| RealVector b1 = previousSol[0].mapMultiply(sigma).subtract(q); | ||
| RealVector b2 = previousSol[2].subtract(Rinv.operate(previousSol[1])); | ||
|
|
||
| //SOLVE KKT SYSYEM | ||
| ADMMQPSolution sol = solve(b1, b2); | ||
| RealVector xtilde = sol.getX(); | ||
| RealVector vtilde = sol.getV(); | ||
|
|
||
| //UPDATE ZTILDE | ||
| RealVector ztilde = zold.add(Rinv.operate(vtilde.subtract(yold))); | ||
| //UPDATE X | ||
| previousSol[0] = xtilde.mapMultiply(alpha).add((xold.mapMultiply(onealfa))); | ||
|
|
||
| //UPDATE Z PARTIAL | ||
| RealVector zpartial = ztilde.mapMultiply(alpha).add(zold.mapMultiply(onealfa)).add(Rinv.operate(yold)); | ||
|
|
||
| //PROJECT ZPARTIAL AND UPDATE Z | ||
| for (int j = 0; j < previousSol[2].getDimension(); j++) { | ||
| previousSol[2].setEntry(j, FastMath.min(FastMath.max(zpartial.getEntry(j), lb.getEntry(j)), ub.getEntry(j))); | ||
| } | ||
|
|
||
| //UPDATE Y | ||
| RealVector ytilde = ztilde.mapMultiply(alpha).add(zold.mapMultiply(onealfa).subtract(previousSol[2])); | ||
| previousSol[1] = yold.add(R.operate(ytilde)); | ||
|
|
||
| return new ADMMQPSolution(previousSol[0], vtilde, previousSol[1], previousSol[2]); | ||
| } | ||
|
|
||
| } |
Oops, something went wrong.