28055112333

Committed 23 Jun 2026 08:32PM UTC coverage: 74.885% (+0.04%) from 74.843%

Build # 28055112333

Build Type

Pull #2635

github

Committed by

web-flow

Commit Message

Merge 777a46dd1 into 9863b578a

Pull Request Pull Request #2635: Add L-BFGS-B Limited-Memory Bound-Constrained Optimizer

Coverage Stats

238 of 285 new or added lines in 1 file covered. (83.51%)

3 existing lines in 2 files now uncovered.

59504 of 79461 relevant lines covered (74.88%)

8626878.93 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

83.51

/ql/math/optimization/lbfgsb.cpp

/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */

/*
 Copyright (C) 2026 Colin Alberts

 This file is part of QuantLib, a free-software/open-source library
 for financial quantitative analysts and developers - http://quantlib.org/

 QuantLib is free software: you can redistribute it and/or modify it
 under the terms of the QuantLib license.  You should have received a
 copy of the license along with this program; if not, please email
 <quantlib-dev@lists.sf.net>. The license is also available online at
 <https://www.quantlib.org/license.shtml>.

 This program is distributed in the hope that it will be useful, but WITHOUT
 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 FOR A PARTICULAR PURPOSE.  See the license for more details.
*/

#include <ql/math/matrix.hpp>
#include <ql/math/optimization/constraint.hpp>
#include <ql/math/optimization/lbfgsb.hpp>
#include <ql/math/optimization/problem.hpp>
#include <algorithm>
#include <cmath>
#include <deque>
#include <utility>
#include <vector>

namespace QuantLib {

    namespace {

        const Real INF = QL_MAX_REAL;

        // A bound is "absent" when it equals +/- DBL_MAX (the value
        // returned by the default Constraint and by NoConstraint).
        // The 0.5 factor guards against overflow in subsequent arithmetic
        // such as (x - bound) or (bound - x).
        inline bool noUpper(Real u) {
            return u >= 0.5 * INF;
        }
        inline bool noLower(Real l) {
            return l <= -0.5 * INF;
        }

        // Compact limited-memory representation of the BFGS Hessian
        // approximation  B = theta I - W M W^T   (Byrd, Lu, Nocedal & Zhu
        // 1995, eq. 3.5), built from the stored correction pairs.
        struct CompactRep {
            Matrix W; // n x 2col,   [ Y | theta S ]
            Matrix M; // 2col x 2col, inverse of the middle matrix
            Real theta = 1.0;
            Size col = 0;
        };

        CompactRep
        buildCompactRep(const std::deque<Array>& S, const std::deque<Array>& Y, Real theta) {
            CompactRep rep;

            rep.theta = theta;
            rep.col = S.size();
            const Size col = rep.col;

            if (col == 0)
                return rep;

            const Size n = S.front().size();
            Matrix W(n, 2 * col);

            for (Size j = 0; j < col; ++j)
                for (Size i = 0; i < n; ++i) {
                    W[i][j] = Y[j][i];
                    W[i][col + j] = theta * S[j][i];
                }

            // small Gram matrices S^T S and S^T Y
            Matrix StS(col, col), StY(col, col);

            for (Size i = 0; i < col; ++i)
                for (Size j = 0; j < col; ++j) {
                    StS[i][j] = DotProduct(S[i], S[j]);
                    StY[i][j] = DotProduct(S[i], Y[j]);
                }

            // middle matrix  [ -D  L^T ; L  theta S^T S ]  with
            // D = diag(s_i.y_i) and L the strictly lower part of S^T Y.
            Matrix mid(2 * col, 2 * col, 0.0);

            for (Size i = 0; i < col; ++i) {
                mid[i][i] = -StY[i][i]; // -D

                for (Size j = 0; j < col; ++j) {
                    if (i > j)
                        mid[col + i][j] = StY[i][j]; // L
                    if (j > i)
                        mid[i][col + j] = StY[j][i]; // L^T

                    mid[col + i][col + j] = theta * StS[i][j];
                }
            }

            rep.W = W;
            rep.M = inverse(mid);
            return rep;
        }

        inline Array wRow(const Matrix& W, Size i) {
            Array r(W.columns());

            for (Size k = 0; k < r.size(); ++k)
                r[k] = W[i][k];

            return r;
        }

        // Generalized Cauchy point (Byrd et al. 1995, Algorithm CP).
        // Walks the breakpoints of the projected steepest-descent path
        // x(t) = P(x - t g, l, u) and locates the first local minimizer of
        // the quadratic model along it.  Returns the Cauchy point xcp, the
        // accumulated c = W^T (xcp - x) needed by the subspace step, and
        // the set of variables left free (not pinned at a bound).
        void generalizedCauchyPoint(const Array& x,
                                    const Array& g,
                                    const Array& lo,
                                    const Array& hi,
                                    const CompactRep& rep,
                                    Array& xcp,
                                    Array& cOut,
                                    std::vector<bool>& isFree) {
            const Size n = x.size();
            const Size m2 = 2 * rep.col;

            xcp = x;
            isFree.assign(n, true);

            Array t(n), d(n, 0.0);
            std::vector<Size> brk; // indices with a strictly positive breakpoint
            for (Size i = 0; i < n; ++i) {
                Real ti;

                if (g[i] < 0.0)
                    ti = noUpper(hi[i]) ? INF : (x[i] - hi[i]) / g[i];
                else if (g[i] > 0.0)
                    ti = noLower(lo[i]) ? INF : (x[i] - lo[i]) / g[i];
                else
                    ti = INF;

                t[i] = ti;

                if (ti <= 0.0) {
                    // already sitting on the bound the gradient pushes into
                    d[i] = 0.0;
                    isFree[i] = false;
                } else {
                    d[i] = -g[i];
                    brk.push_back(i);
                }
            }

            Array c(m2, 0.0);
            if (brk.empty()) { // every variable is pinned
                cOut = c;
                return;
            }

            std::sort(brk.begin(), brk.end(), [&t](Size a, Size b) { return t[a] < t[b]; });

            Array p = (rep.col > 0) ? (d * rep.W) : Array(0);

            Real fp = -DotProduct(d, d); // first derivative of the model, m'(0)
            Real fpp = -rep.theta * fp;  // second derivative, theta d^T d - ...

            if (rep.col > 0)
                fpp -= DotProduct(p, rep.M * p);
            Real fppFloor = QL_EPSILON * (fpp > 0.0 ? fpp : 1.0);

            Real dtMin = (fpp > 0.0) ? -fp / fpp : 0.0;

            Real tOld = 0.0;
            Size ptr = 0;
            Size b = brk[ptr];
            Real tb = t[b];
            Real dt = tb;
            bool exhausted = false;

            while (dtMin >= dt && tb < INF) {
                // pin variable b at the bound it has reached
                Real xcpb = (d[b] > 0.0) ? hi[b] : (d[b] < 0.0 ? lo[b] : x[b]);
                xcp[b] = xcpb;
                isFree[b] = false;
                Real zb = xcpb - x[b];
                Real gb = g[b];

                if (rep.col > 0) {
                    c += dt * p;

                    Array wb = wRow(rep.W, b);
                    Array Mc = rep.M * c;
                    Array Mp = rep.M * p;
                    Array Mwb = rep.M * wb;

                    fp += dt * fpp + gb * gb + rep.theta * gb * zb - gb * DotProduct(wb, Mc);
                    fpp += -rep.theta * gb * gb - 2.0 * gb * DotProduct(wb, Mp) -
                           gb * gb * DotProduct(wb, Mwb);
                    p += gb * wb;
                } else {
                    fp += dt * fpp + gb * gb + rep.theta * gb * zb;
                    fpp += -rep.theta * gb * gb;
                }

                if (fpp < fppFloor)
                    fpp = fppFloor;

                d[b] = 0.0;
                dtMin = (fpp > 0.0) ? -fp / fpp : 0.0;
                tOld = tb;

                if (++ptr >= brk.size()) {
                    exhausted = true;
                    break;
                }

                b = brk[ptr];
                tb = t[b];
                dt = tb - tOld;
            }

            // advance the still-free variables along the final segment
            if (exhausted)
                dtMin = 0.0;

            dtMin = std::max(dtMin, Real(0.0));
            tOld += dtMin;

            for (Size i = 0; i < n; ++i)
                if (isFree[i])
                    xcp[i] = x[i] + tOld * d[i];

            if (rep.col > 0)
                c += dtMin * p;
            cOut = c;
        }

        // Subspace minimization by the direct primal method (Byrd et al.
        // 1995, section 5.1).  Minimizes the quadratic model over the free
        // variables, holding the active ones at their Cauchy-point bound,
        // then truncates the step back into the box.
        Array subspaceMinimization(const Array& x,
                                   const Array& g,
                                   const Array& xcp,
                                   const Array& c,
                                   const Array& lo,
                                   const Array& hi,
                                   const CompactRep& rep,
                                   const std::vector<bool>& isFree) {
            const Size n = x.size();
            std::vector<Size> freeIdx;

            for (Size i = 0; i < n; ++i)
                if (isFree[i])
                    freeIdx.push_back(i);

            const Size nf = freeIdx.size();

            Array xbar = xcp;
            if (nf == 0)
                return xbar;

            const Real theta = rep.theta;
            const Size m2 = 2 * rep.col;

            // reduced gradient of the model at the Cauchy point:
            //   r = [ g + theta (xcp - x) - W M c ]  restricted to free vars
            Array WMc(n, 0.0);
            if (rep.col > 0)
                WMc = rep.W * (rep.M * c);

            Array r(nf);
            for (Size k = 0; k < nf; ++k) {
                Size i = freeIdx[k];
                r[k] = g[i] + theta * (xcp[i] - x[i]) - WMc[i];
            }

            Array dhat(nf);
            if (rep.col == 0) {
                for (Size k = 0; k < nf; ++k)
                    dhat[k] = -r[k] / theta;
            } else {
                // v = M (W_free^T r)
                Array wtr(m2, 0.0);

                for (Size k = 0; k < nf; ++k) {
                    Size i = freeIdx[k];

                    for (Size j = 0; j < m2; ++j)
                        wtr[j] += rep.W[i][j] * r[k];
                }

                Array v = rep.M * wtr;

                // N = I - (1/theta) M (W_free^T W_free)
                Matrix WftWf(m2, m2, 0.0);

                for (Size k = 0; k < nf; ++k) {
                    Size i = freeIdx[k];

                    for (Size a = 0; a < m2; ++a)
                        for (Size bb = 0; bb < m2; ++bb)
                            WftWf[a][bb] += rep.W[i][a] * rep.W[i][bb];
                }

                Matrix N = rep.M * WftWf;
                N *= -1.0 / theta;

                for (Size a = 0; a < m2; ++a)
                    N[a][a] += 1.0;
                v = inverse(N) * v;

                // dhat = -(1/theta) r - (1/theta^2) W_free v
                for (Size k = 0; k < nf; ++k) {
                    Size i = freeIdx[k];
                    Real Wfv = 0.0;

                    for (Size j = 0; j < m2; ++j)
                        Wfv += rep.W[i][j] * v[j];

                    dhat[k] = -r[k] / theta - Wfv / (theta * theta);
                }
            }

            // truncate so that every free variable stays inside the box
            Real alphaStar = 1.0;

            for (Size k = 0; k < nf; ++k) {
                Size i = freeIdx[k];

                if (dhat[k] > 0.0 && !noUpper(hi[i]))
                    alphaStar = std::min(alphaStar, (hi[i] - xcp[i]) / dhat[k]);
                else if (dhat[k] < 0.0 && !noLower(lo[i]))
                    alphaStar = std::min(alphaStar, (lo[i] - xcp[i]) / dhat[k]);
            }
            alphaStar = std::max(alphaStar, Real(0.0));

            for (Size k = 0; k < nf; ++k)
                xbar[freeIdx[k]] = xcp[freeIdx[k]] + alphaStar * dhat[k];
            return xbar;
        }

        // Largest feasible step length along d starting from x.
        Real maxFeasibleStep(const Array& x, const Array& d, const Array& lo, const Array& hi) {
            Real stp = INF;

            for (Size i = 0; i < x.size(); ++i) {
                if (d[i] > 0.0 && !noUpper(hi[i]))
                    stp = std::min(stp, (hi[i] - x[i]) / d[i]);
                else if (d[i] < 0.0 && !noLower(lo[i]))
                    stp = std::min(stp, (lo[i] - x[i]) / d[i]);
            }
            return stp;
        }

        // Line search enforcing the strong Wolfe conditions (Nocedal &
        // Wright, Algorithms 3.5/3.6), capped at the largest feasible step.
        // On success returns the step, trial point, value and gradient.
        bool lineSearchWolfe(Problem& P,
                             const Array& x,
                             const Array& d,
                             Real f0,
                             const Array& g0,
                             Real stpMax,
                             Real& alpha,
                             Array& xt,
                             Real& ft,
                             Array& gt) {
            const Real c1 = 1e-4, c2 = 0.9;
            const Size maxIter = 30;

            Real dphi0 = DotProduct(g0, d);
            if (dphi0 >= 0.0)
                return false; // not a descent direction

            bool haveBest = false;
            Real bestAlpha = 0.0, bestF = f0;
            Array bestX, bestG;

            // evaluate phi(a) and phi'(a), tracking the best decrease seen
            auto eval = [&](Real a) -> Real {
                xt = x + a * d;
                ft = P.valueAndGradient(gt, xt);

                if (ft < bestF) {
                    haveBest = true;
                    bestF = ft;
                    bestAlpha = a;
                    bestX = xt;
                    bestG = gt;
                }
                return DotProduct(gt, d);
            };

            Real aLo = 0.0, aHi = 0.0, fLo = f0;
            bool bracketed = false;

            Real aPrev = 0.0, fPrev = f0;
            Real a = std::min(Real(1.0), stpMax);

            for (Size i = 0; i < maxIter; ++i) {
                Real dphi = eval(a);

                if (ft > f0 + c1 * a * dphi0 || (i > 0 && ft >= fPrev)) {
                    aLo = aPrev;
                    fLo = fPrev;
                    aHi = a;
                    bracketed = true;
                    break;
                }

                if (std::fabs(dphi) <= -c2 * dphi0) {
                    alpha = a;
                    return true; // strong Wolfe satisfied
                }

                if (dphi >= 0.0) {
                    aLo = a;
                    fLo = ft;
                    aHi = aPrev;
                    bracketed = true;
                    break;
                }

                aPrev = a;
                fPrev = ft;

                if (a >= stpMax)
                    break; // cannot expand further
                a = std::min(Real(2.0) * a, stpMax);
            }

            if (bracketed) {
                for (Size j = 0; j < maxIter; ++j) {
                    a = 0.5 * (aLo + aHi); // bisection
                    Real dphi = eval(a);

                    if (ft > f0 + c1 * a * dphi0 || ft >= fLo) {
                        aHi = a;
                    } else {
                        if (std::fabs(dphi) <= -c2 * dphi0) {
                            alpha = a;
                            return true;
                        }

                        if (dphi * (aHi - aLo) >= 0.0)
                            aHi = aLo;

                        aLo = a;
                        fLo = ft;
                    }

                    if (std::fabs(aHi - aLo) < QL_EPSILON * std::max(Real(1.0), std::fabs(a)))
                        break;
                }
            }

            // strong Wolfe not reached: accept the best sufficient decrease
            if (haveBest) {
                alpha = bestAlpha;
                xt = bestX;
                ft = bestF;
                gt = bestG;

                return true;
            }
            return false;
        }

    }

    LBFGSB::LBFGSB(Size memory, Real pgTol, Real fTol)
    : m_(memory), pgTol_(pgTol), factr_(fTol / QL_EPSILON) {
        QL_REQUIRE(memory > 0, "memory must be positive");
    }

    LBFGSB::LBFGSB(Array lowerBound, Array upperBound, Size memory, Real pgTol, Real fTol)
    : m_(memory), pgTol_(pgTol), factr_(fTol / QL_EPSILON), lowerBound_(std::move(lowerBound)),
      upperBound_(std::move(upperBound)) {
        QL_REQUIRE(memory > 0, "memory must be positive");
        QL_REQUIRE(lowerBound_.size() == upperBound_.size(),
                   "lower and upper bound sizes are inconsistent");
    }

    EndCriteria::Type LBFGSB::minimize(Problem& P, const EndCriteria& endCriteria) {
        EndCriteria::Type ecType = EndCriteria::None;
        P.reset();

        Array x = P.currentValue();
        const Size n = x.size();

        Array lo = lowerBound_.empty() ? P.constraint().lowerBound(x) : lowerBound_;
        Array hi = upperBound_.empty() ? P.constraint().upperBound(x) : upperBound_;

        QL_REQUIRE(lo.size() == n && hi.size() == n,
                   "bounds size does not match the number of variables");

        // start from a feasible point
        for (Size i = 0; i < n; ++i)
            x[i] = std::min(std::max(x[i], lo[i]), hi[i]);

        Array g(n);
        Real f = P.valueAndGradient(g, x);
        P.setCurrentValue(x);
        P.setFunctionValue(f);

        std::deque<Array> S, Y;
        Real theta = 1.0;
        Size iter = 0;
        Real pgInf = 0.0; // projected gradient ∞-norm; persists for final reporting

        while (true) {
            // infinity norm of the projected gradient
            pgInf = 0.0;

            for (Size i = 0; i < n; ++i) {
                Real proj = std::min(std::max(x[i] - g[i], lo[i]), hi[i]) - x[i];
                pgInf = std::max(pgInf, std::fabs(proj));
            }

            P.setGradientNormValue(pgInf * pgInf);

            if (pgInf < pgTol_) {
                ecType = EndCriteria::ZeroGradientNorm;
                break;
            }

            if (endCriteria.checkZeroGradientNorm(pgInf, ecType))
                break;

            if (endCriteria.checkMaxIterations(iter, ecType))
                break;

            // compact representation; drop the oldest pairs if the middle
            // matrix turns out numerically singular
            CompactRep rep;

            while (true) {
                try {
                    rep = buildCompactRep(S, Y, theta);
                    break;
                } catch (...) {
                    if (S.empty()) {
                        rep = CompactRep();
                        rep.theta = theta;
                        break;
                    }
                    S.pop_front();
                    Y.pop_front();
                }
            }

            Array xcp, c;
            std::vector<bool> isFree;
            generalizedCauchyPoint(x, g, lo, hi, rep, xcp, c, isFree);

            Array xbar;
            try {
                xbar = subspaceMinimization(x, g, xcp, c, lo, hi, rep, isFree);
            } catch (...) {
                xbar = xcp; // fall back to the Cauchy point
            }

            Array d = xbar - x;
            Real dphi0 = DotProduct(g, d);

            // fall back to the (always-descent) Cauchy direction, then to
            // the projected gradient, if the subspace step is not downhill
            if (Norm2(d) < QL_EPSILON || dphi0 >= 0.0) {
                d = xcp - x;
                dphi0 = DotProduct(g, d);
            }

            if (Norm2(d) < QL_EPSILON) {
                ecType = EndCriteria::StationaryPoint;
                break;
            }

            if (dphi0 >= 0.0) {
                for (Size i = 0; i < n; ++i)
                    d[i] = std::min(std::max(x[i] - g[i], lo[i]), hi[i]) - x[i];
                dphi0 = DotProduct(g, d);
                if (dphi0 >= 0.0 || Norm2(d) < QL_EPSILON) {
                    ecType = EndCriteria::StationaryPoint;
                    break;
                }
            }

            Real stpMax = maxFeasibleStep(x, d, lo, hi);
            if (stpMax < QL_EPSILON) {
                ecType = EndCriteria::StationaryPoint;
                break;
            }

            Real alpha = 0.0;
            Array xt(n), gt(n);
            Real ft = 0.0;
            if (!lineSearchWolfe(P, x, d, f, g, stpMax, alpha, xt, ft, gt)) {
                ecType = EndCriteria::StationaryFunctionValue;
                break;
            }

            // limited-memory update with the curvature safeguard
            Array s = xt - x;
            Array y = gt - g;
            Real sy = DotProduct(s, y);
            Real yy = DotProduct(y, y);

            if (yy > 0.0 && sy > QL_EPSILON * yy) {
                S.push_back(s);
                Y.push_back(y);

                if (S.size() > m_) {
                    S.pop_front();
                    Y.pop_front();
                }
                theta = yy / sy;
            }

            Real fOld = f;
            x = xt;
            f = ft;
            g = gt;
            ++iter;

            P.setCurrentValue(x);
            P.setFunctionValue(f);

            // relative function-reduction stop (SciPy's factr criterion)
            Real denom = std::max(std::max(std::fabs(fOld), std::fabs(f)), Real(1.0));
            if ((fOld - f) <= factr_ * QL_EPSILON * denom) {
                ecType = EndCriteria::StationaryFunctionValue;
                break;
            }
        }

        P.setCurrentValue(x);
        P.setFunctionValue(f);
        P.setGradientNormValue(pgInf * pgInf);
        return ecType;
    }

}

1	/* -- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -- */
2
3	/*
4	Copyright (C) 2026 Colin Alberts
5
6	This file is part of QuantLib, a free-software/open-source library
7	for financial quantitative analysts and developers - http://quantlib.org/
8
9	QuantLib is free software: you can redistribute it and/or modify it
10	under the terms of the QuantLib license. You should have received a
11	copy of the license along with this program; if not, please email
12	<quantlib-dev@lists.sf.net>. The license is also available online at
13	<https://www.quantlib.org/license.shtml>.
14
15	This program is distributed in the hope that it will be useful, but WITHOUT
16	ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17	FOR A PARTICULAR PURPOSE. See the license for more details.
18	*/
19
20	#include <ql/math/matrix.hpp>
21	#include <ql/math/optimization/constraint.hpp>
22	#include <ql/math/optimization/lbfgsb.hpp>
23	#include <ql/math/optimization/problem.hpp>
24	#include <algorithm>
25	#include <cmath>
26	#include <deque>
27	#include <utility>
28	#include <vector>
29
30	namespace QuantLib {
31
32	namespace {
33
34	const Real INF = QL_MAX_REAL;
35
36	// A bound is "absent" when it equals +/- DBL_MAX (the value
37	// returned by the default Constraint and by NoConstraint).
38	// The 0.5 factor guards against overflow in subsequent arithmetic
39	// such as (x - bound) or (bound - x).
40	inline bool noUpper(Real u) {
41	return u >= 0.5 * INF;
42	}
43	inline bool noLower(Real l) {
44	return l <= -0.5 * INF;
45	}
46
47	// Compact limited-memory representation of the BFGS Hessian
48	// approximation B = theta I - W M W^T (Byrd, Lu, Nocedal & Zhu
49	// 1995, eq. 3.5), built from the stored correction pairs.
50	struct CompactRep {	193✔
51	Matrix W; // n x 2col, [ Y \| theta S ]
52	Matrix M; // 2col x 2col, inverse of the middle matrix
53	Real theta = 1.0;
54	Size col = 0;
55	};
56
57	CompactRep
58	buildCompactRep(const std::deque<Array>& S, const std::deque<Array>& Y, Real theta) {	193✔
59	CompactRep rep;
60
61	rep.theta = theta;	193✔
62	rep.col = S.size();	193✔
63	const Size col = rep.col;
64
65	if (col == 0)	193✔
66	return rep;
67
68	const Size n = S.front().size();
69	Matrix W(n, 2 * col);	180✔
70
71	for (Size j = 0; j < col; ++j)	1,129✔
72	for (Size i = 0; i < n; ++i) {	10,597✔
73	W[i][j] = Y[j][i];	9,648✔
74	W[i][col + j] = theta * S[j][i];	9,648✔
75	}
76
77	// small Gram matrices S^T S and S^T Y
78	Matrix StS(col, col), StY(col, col);	180✔
79
80	for (Size i = 0; i < col; ++i)	1,129✔
81	for (Size j = 0; j < col; ++j) {	7,894✔
82	StS[i][j] = DotProduct(S[i], S[j]);	6,945✔
83	StY[i][j] = DotProduct(S[i], Y[j]);	6,945✔
84	}
85
86	// middle matrix [ -D L^T ; L theta S^T S ] with
87	// D = diag(s_i.y_i) and L the strictly lower part of S^T Y.
88	Matrix mid(2 * col, 2 * col, 0.0);	180✔
89
90	for (Size i = 0; i < col; ++i) {	1,129✔
91	mid[i][i] = -StY[i][i]; // -D	949✔
92
93	for (Size j = 0; j < col; ++j) {	7,894✔
94	if (i > j)	6,945✔
95	mid[col + i][j] = StY[i][j]; // L	2,998✔
96	if (j > i)	6,945✔
97	mid[i][col + j] = StY[j][i]; // L^T	2,998✔
98
99	mid[col + i][col + j] = theta * StS[i][j];	6,945✔
100	}
101	}
102
103	rep.W = W;	180✔
104	rep.M = inverse(mid);	360✔
105	return rep;
NEW 106	}	×
107
NEW 108	inline Array wRow(const Matrix& W, Size i) {	×
NEW 109	Array r(W.columns());	×
110
NEW 111	for (Size k = 0; k < r.size(); ++k)	×
NEW 112	r[k] = W[i][k];	×
113
NEW 114	return r;	×
115	}
116
117	// Generalized Cauchy point (Byrd et al. 1995, Algorithm CP).
118	// Walks the breakpoints of the projected steepest-descent path
119	// x(t) = P(x - t g, l, u) and locates the first local minimizer of
120	// the quadratic model along it. Returns the Cauchy point xcp, the
121	// accumulated c = W^T (xcp - x) needed by the subspace step, and
122	// the set of variables left free (not pinned at a bound).
123	void generalizedCauchyPoint(const Array& x,	193✔
124	const Array& g,
125	const Array& lo,
126	const Array& hi,
127	const CompactRep& rep,
128	Array& xcp,
129	Array& cOut,
130	std::vector<bool>& isFree) {
131	const Size n = x.size();
132	const Size m2 = 2 * rep.col;	193✔
133
134	xcp = x;	193✔
135	isFree.assign(n, true);
136
137	Array t(n), d(n, 0.0);	193✔
138	std::vector<Size> brk; // indices with a strictly positive breakpoint
139	for (Size i = 0; i < n; ++i) {	2,533✔
140	Real ti;
141
142	if (g[i] < 0.0)	2,340✔
143	ti = noUpper(hi[i]) ? INF : (x[i] - hi[i]) / g[i];	1,162✔
144	else if (g[i] > 0.0)	1,178✔
145	ti = noLower(lo[i]) ? INF : (x[i] - lo[i]) / g[i];	1,177✔
146	else
147	ti = INF;
148
149	t[i] = ti;	2,340✔
150
151	if (ti <= 0.0) {	2,340✔
152	// already sitting on the bound the gradient pushes into
153	d[i] = 0.0;	20✔
154	isFree[i] = false;
155	} else {
156	d[i] = -g[i];	2,320✔
157	brk.push_back(i);	2,320✔
158	}
159	}
160
161	Array c(m2, 0.0);	193✔
162	if (brk.empty()) { // every variable is pinned	193✔
NEW 163	cOut = c;	×
164	return;
165	}
166
167	std::sort(brk.begin(), brk.end(), [&t](Size a, Size b) { return t[a] < t[b]; });	5,733✔
168
169	Array p = (rep.col > 0) ? (d * rep.W) : Array(0);	193✔
170
171	Real fp = -DotProduct(d, d); // first derivative of the model, m'(0)	193✔
172	Real fpp = -rep.theta * fp; // second derivative, theta d^T d - ...	193✔
173
174	if (rep.col > 0)	193✔
175	fpp -= DotProduct(p, rep.M * p);	360✔
176	Real fppFloor = QL_EPSILON * (fpp > 0.0 ? fpp : 1.0);	193✔
177
178	Real dtMin = (fpp > 0.0) ? -fp / fpp : 0.0;	193✔
179
180	Real tOld = 0.0;
181	Size ptr = 0;
182	Size b = brk[ptr];	193✔
183	Real tb = t[b];	193✔
184	Real dt = tb;
185	bool exhausted = false;
186
187	while (dtMin >= dt && tb < INF) {	204✔
188	// pin variable b at the bound it has reached
189	Real xcpb = (d[b] > 0.0) ? hi[b] : (d[b] < 0.0 ? lo[b] : x[b]);	15✔
190	xcp[b] = xcpb;	15✔
191	isFree[b] = false;
192	Real zb = xcpb - x[b];	15✔
193	Real gb = g[b];
194
195	if (rep.col > 0) {	15✔
NEW 196	c += dt * p;	×
197
NEW 198	Array wb = wRow(rep.W, b);	×
NEW 199	Array Mc = rep.M * c;	×
NEW 200	Array Mp = rep.M * p;	×
NEW 201	Array Mwb = rep.M * wb;	×
202
NEW 203	fp += dt * fpp + gb * gb + rep.theta * gb * zb - gb * DotProduct(wb, Mc);	×
NEW 204	fpp += -rep.theta * gb * gb - 2.0 * gb * DotProduct(wb, Mp) -	×
NEW 205	gb * gb * DotProduct(wb, Mwb);	×
NEW 206	p += gb * wb;	×
207	} else {
208	fp += dt * fpp + gb * gb + rep.theta * gb * zb;	15✔
209	fpp += -rep.theta * gb * gb;	15✔
210	}
211
212	if (fpp < fppFloor)	15✔
213	fpp = fppFloor;
214
215	d[b] = 0.0;	15✔
216	dtMin = (fpp > 0.0) ? -fp / fpp : 0.0;	15✔
217	tOld = tb;
218
219	if (++ptr >= brk.size()) {	15✔
220	exhausted = true;
221	break;
222	}
223
224	b = brk[ptr];	11✔
225	tb = t[b];	11✔
226	dt = tb - tOld;	11✔
227	}
228
229	// advance the still-free variables along the final segment
230	if (exhausted)	193✔
231	dtMin = 0.0;	4✔
232
233	dtMin = std::max(dtMin, Real(0.0));	193✔
234	tOld += dtMin;	193✔
235
236	for (Size i = 0; i < n; ++i)	2,533✔
237	if (isFree[i])	2,340✔
238	xcp[i] = x[i] + tOld * d[i];	2,305✔
239
240	if (rep.col > 0)	193✔
241	c += dtMin * p;	360✔
242	cOut = c;	193✔
243	}
244
245	// Subspace minimization by the direct primal method (Byrd et al.
246	// 1995, section 5.1). Minimizes the quadratic model over the free
247	// variables, holding the active ones at their Cauchy-point bound,
248	// then truncates the step back into the box.
249	Array subspaceMinimization(const Array& x,	193✔
250	const Array& g,
251	const Array& xcp,
252	const Array& c,
253	const Array& lo,
254	const Array& hi,
255	const CompactRep& rep,
256	const std::vector<bool>& isFree) {
257	const Size n = x.size();
258	std::vector<Size> freeIdx;
259
260	for (Size i = 0; i < n; ++i)	2,533✔
261	if (isFree[i])	2,340✔
262	freeIdx.push_back(i);	2,305✔
263
264	const Size nf = freeIdx.size();
265
266	Array xbar = xcp;	193✔
267	if (nf == 0)	193✔
268	return xbar;
269
270	const Real theta = rep.theta;	189✔
271	const Size m2 = 2 * rep.col;	189✔
272
273	// reduced gradient of the model at the Cauchy point:
274	// r = [ g + theta (xcp - x) - W M c ] restricted to free vars
275	Array WMc(n, 0.0);	189✔
276	if (rep.col > 0)	189✔
277	WMc = rep.W * (rep.M * c);	540✔
278
279	Array r(nf);	189✔
280	for (Size k = 0; k < nf; ++k) {	2,494✔
281	Size i = freeIdx[k];	2,305✔
282	r[k] = g[i] + theta * (xcp[i] - x[i]) - WMc[i];	2,305✔
283	}
284
285	Array dhat(nf);	189✔
286	if (rep.col == 0) {	189✔
287	for (Size k = 0; k < nf; ++k)	48✔
288	dhat[k] = -r[k] / theta;	39✔
289	} else {
290	// v = M (W_free^T r)
291	Array wtr(m2, 0.0);	180✔
292
293	for (Size k = 0; k < nf; ++k) {	2,446✔
294	Size i = freeIdx[k];	2,266✔
295
296	for (Size j = 0; j < m2; ++j)	21,460✔
297	wtr[j] += rep.W[i][j] * r[k];	19,194✔
298	}
299
300	Array v = rep.M * wtr;	180✔
301
302	// N = I - (1/theta) M (W_free^T W_free)
303	Matrix WftWf(m2, m2, 0.0);	180✔
304
305	for (Size k = 0; k < nf; ++k) {	2,446✔
306	Size i = freeIdx[k];	2,266✔
307
308	for (Size a = 0; a < m2; ++a)	21,460✔
309	for (Size bb = 0; bb < m2; ++bb)	246,790✔
310	WftWf[a][bb] += rep.W[i][a] * rep.W[i][bb];	227,596✔
311	}
312
313	Matrix N = rep.M * WftWf;	180✔
314	N *= -1.0 / theta;	180✔
315
316	for (Size a = 0; a < m2; ++a)	2,078✔
317	N[a][a] += 1.0;	1,898✔
318	v = inverse(N) * v;	360✔
319
320	// dhat = -(1/theta) r - (1/theta^2) W_free v
321	for (Size k = 0; k < nf; ++k) {	2,446✔
322	Size i = freeIdx[k];	2,266✔
323	Real Wfv = 0.0;
324
325	for (Size j = 0; j < m2; ++j)	21,460✔
326	Wfv += rep.W[i][j] * v[j];	19,194✔
327
328	dhat[k] = -r[k] / theta - Wfv / (theta * theta);	2,266✔
329	}
330	}
331
332	// truncate so that every free variable stays inside the box
333	Real alphaStar = 1.0;	189✔
334
335	for (Size k = 0; k < nf; ++k) {	2,494✔
336	Size i = freeIdx[k];	2,305✔
337
338	if (dhat[k] > 0.0 && !noUpper(hi[i]))	2,305✔
339	alphaStar = std::min(alphaStar, (hi[i] - xcp[i]) / dhat[k]);	29✔
340	else if (dhat[k] < 0.0 && !noLower(lo[i]))	2,276✔
341	alphaStar = std::min(alphaStar, (lo[i] - xcp[i]) / dhat[k]);	26✔
342	}
343	alphaStar = std::max(alphaStar, Real(0.0));	189✔
344
345	for (Size k = 0; k < nf; ++k)	2,494✔
346	xbar[freeIdx[k]] = xcp[freeIdx[k]] + alphaStar * dhat[k];	2,305✔
347	return xbar;
348	}
349
350	// Largest feasible step length along d starting from x.
351	Real maxFeasibleStep(const Array& x, const Array& d, const Array& lo, const Array& hi) {	193✔
352	Real stp = INF;	193✔
353
354	for (Size i = 0; i < x.size(); ++i) {	2,533✔
355	if (d[i] > 0.0 && !noUpper(hi[i]))	2,340✔
356	stp = std::min(stp, (hi[i] - x[i]) / d[i]);	71✔
357	else if (d[i] < 0.0 && !noLower(lo[i]))	2,295✔
358	stp = std::min(stp, (lo[i] - x[i]) / d[i]);	49✔
359	}
360	return stp;	193✔
361	}
362
363	// Line search enforcing the strong Wolfe conditions (Nocedal &
364	// Wright, Algorithms 3.5/3.6), capped at the largest feasible step.
365	// On success returns the step, trial point, value and gradient.
366	bool lineSearchWolfe(Problem& P,	193✔
367	const Array& x,
368	const Array& d,
369	Real f0,
370	const Array& g0,
371	Real stpMax,
372	Real& alpha,
373	Array& xt,
374	Real& ft,
375	Array& gt) {
376	const Real c1 = 1e-4, c2 = 0.9;
377	const Size maxIter = 30;
378
379	Real dphi0 = DotProduct(g0, d);	193✔
380	if (dphi0 >= 0.0)	193✔
381	return false; // not a descent direction
382
383	bool haveBest = false;	193✔
384	Real bestAlpha = 0.0, bestF = f0;	193✔
385	Array bestX, bestG;
386
387	// evaluate phi(a) and phi'(a), tracking the best decrease seen
388	auto eval = [&](Real a) -> Real {	246✔
389	xt = x + a * d;	492✔
390	ft = P.valueAndGradient(gt, xt);	246✔
391
392	if (ft < bestF) {	246✔
393	haveBest = true;	199✔
394	bestF = ft;	199✔
395	bestAlpha = a;	199✔
396	bestX = xt;	199✔
397	bestG = gt;	199✔
398	}
399	return DotProduct(gt, d);	246✔
400	};	193✔
401
402	Real aLo = 0.0, aHi = 0.0, fLo = f0;
403	bool bracketed = false;
404
405	Real aPrev = 0.0, fPrev = f0;
406	Real a = std::min(Real(1.0), stpMax);	193✔
407
408	for (Size i = 0; i < maxIter; ++i) {	197✔
409	Real dphi = eval(a);	197✔
410
411	if (ft > f0 + c1 * a * dphi0 \|\| (i > 0 && ft >= fPrev)) {	197✔
412	aLo = aPrev;
413	fLo = fPrev;
414	aHi = a;
415	bracketed = true;
416	break;
417	}
418
419	if (std::fabs(dphi) <= -c2 * dphi0) {	177✔
420	alpha = a;	170✔
421	return true; // strong Wolfe satisfied	170✔
422	}
423
424	if (dphi >= 0.0) {	7✔
425	aLo = a;
426	fLo = ft;
427	aHi = aPrev;
428	bracketed = true;
429	break;
430	}
431
432	aPrev = a;
433	fPrev = ft;
434
435	if (a >= stpMax)	5✔
436	break; // cannot expand further
437	a = std::min(Real(2.0) * a, stpMax);	4✔
438	}
439
440	if (bracketed) {	23✔
441	for (Size j = 0; j < maxIter; ++j) {	49✔
442	a = 0.5 * (aLo + aHi); // bisection	49✔
443	Real dphi = eval(a);	49✔
444
445	if (ft > f0 + c1 * a * dphi0 \|\| ft >= fLo) {	49✔
446	aHi = a;
447	} else {
448	if (std::fabs(dphi) <= -c2 * dphi0) {	22✔
449	alpha = a;	22✔
450	return true;	22✔
451	}
452
NEW 453	if (dphi * (aHi - aLo) >= 0.0)	×
454	aHi = aLo;
455
456	aLo = a;
457	fLo = ft;
458	}
459
460	if (std::fabs(aHi - aLo) < QL_EPSILON * std::max(Real(1.0), std::fabs(a)))	27✔
461	break;
462	}
463	}
464
465	// strong Wolfe not reached: accept the best sufficient decrease
466	if (haveBest) {	1✔
467	alpha = bestAlpha;	1✔
468	xt = bestX;	1✔
469	ft = bestF;	1✔
470	gt = bestG;	1✔
471
472	return true;
473	}
474	return false;
475	}
476
477	}
478
479	LBFGSB::LBFGSB(Size memory, Real pgTol, Real fTol)	13✔
480	: m_(memory), pgTol_(pgTol), factr_(fTol / QL_EPSILON) {	13✔
481	QL_REQUIRE(memory > 0, "memory must be positive");	13✔
482	}	13✔
483
NEW 484	LBFGSB::LBFGSB(Array lowerBound, Array upperBound, Size memory, Real pgTol, Real fTol)	×
NEW 485	: m_(memory), pgTol_(pgTol), factr_(fTol / QL_EPSILON), lowerBound_(std::move(lowerBound)),	×
NEW 486	upperBound_(std::move(upperBound)) {	×
NEW 487	QL_REQUIRE(memory > 0, "memory must be positive");	×
NEW 488	QL_REQUIRE(lowerBound_.size() == upperBound_.size(),	×
489	"lower and upper bound sizes are inconsistent");
NEW 490	}	×
491
492	EndCriteria::Type LBFGSB::minimize(Problem& P, const EndCriteria& endCriteria) {	13✔
493	EndCriteria::Type ecType = EndCriteria::None;	13✔
494	P.reset();
495
496	Array x = P.currentValue();	13✔
497	const Size n = x.size();
498
499	Array lo = lowerBound_.empty() ? P.constraint().lowerBound(x) : lowerBound_;	13✔
500	Array hi = upperBound_.empty() ? P.constraint().upperBound(x) : upperBound_;	13✔
501
502	QL_REQUIRE(lo.size() == n && hi.size() == n,	13✔
503	"bounds size does not match the number of variables");
504
505	// start from a feasible point
506	for (Size i = 0; i < n; ++i)	69✔
507	x[i] = std::min(std::max(x[i], lo[i]), hi[i]);	56✔
508
509	Array g(n);	13✔
510	Real f = P.valueAndGradient(g, x);	13✔
511	P.setCurrentValue(x);	26✔
512	P.setFunctionValue(f);
513
514	std::deque<Array> S, Y;
515	Real theta = 1.0;
516	Size iter = 0;
517	Real pgInf = 0.0; // projected gradient ∞-norm; persists for final reporting
518
519	while (true) {
520	// infinity norm of the projected gradient
521	pgInf = 0.0;	201✔
522
523	for (Size i = 0; i < n; ++i) {	2,559✔
524	Real proj = std::min(std::max(x[i] - g[i], lo[i]), hi[i]) - x[i];	2,358✔
525	pgInf = std::max(pgInf, std::fabs(proj));	2,986✔
526	}
527
528	P.setGradientNormValue(pgInf * pgInf);	201✔
529
530	if (pgInf < pgTol_) {	201✔
531	ecType = EndCriteria::ZeroGradientNorm;	8✔
532	break;	8✔
533	}
534
535	if (endCriteria.checkZeroGradientNorm(pgInf, ecType))	193✔
536	break;
537
538	if (endCriteria.checkMaxIterations(iter, ecType))	193✔
539	break;
540
541	// compact representation; drop the oldest pairs if the middle
542	// matrix turns out numerically singular
543	CompactRep rep;
544
545	while (true) {
546	try {
547	rep = buildCompactRep(S, Y, theta);	193✔
548	break;	193✔
NEW 549	} catch (...) {	×
NEW 550	if (S.empty()) {	×
NEW 551	rep = CompactRep();	×
NEW 552	rep.theta = theta;	×
553	break;
554	}
NEW 555	S.pop_front();	×
NEW 556	Y.pop_front();	×
NEW 557	}	×
558	}
559
560	Array xcp, c;
561	std::vector<bool> isFree;
562	generalizedCauchyPoint(x, g, lo, hi, rep, xcp, c, isFree);	193✔
563
564	Array xbar;
565	try {
566	xbar = subspaceMinimization(x, g, xcp, c, lo, hi, rep, isFree);	386✔
NEW 567	} catch (...) {	×
NEW 568	xbar = xcp; // fall back to the Cauchy point	×
NEW 569	}	×
570
571	Array d = xbar - x;	193✔
572	Real dphi0 = DotProduct(g, d);	193✔
573
574	// fall back to the (always-descent) Cauchy direction, then to
575	// the projected gradient, if the subspace step is not downhill
576	if (Norm2(d) < QL_EPSILON \|\| dphi0 >= 0.0) {	193✔
NEW 577	d = xcp - x;	×
NEW 578	dphi0 = DotProduct(g, d);	×
579	}
580
581	if (Norm2(d) < QL_EPSILON) {	193✔
NEW 582	ecType = EndCriteria::StationaryPoint;	×
NEW 583	break;	×
584	}
585
586	if (dphi0 >= 0.0) {	193✔
NEW 587	for (Size i = 0; i < n; ++i)	×
NEW 588	d[i] = std::min(std::max(x[i] - g[i], lo[i]), hi[i]) - x[i];	×
NEW 589	dphi0 = DotProduct(g, d);	×
NEW 590	if (dphi0 >= 0.0 \|\| Norm2(d) < QL_EPSILON) {	×
NEW 591	ecType = EndCriteria::StationaryPoint;	×
NEW 592	break;	×
593	}
594	}
595
596	Real stpMax = maxFeasibleStep(x, d, lo, hi);	193✔
597	if (stpMax < QL_EPSILON) {	193✔
NEW 598	ecType = EndCriteria::StationaryPoint;	×
NEW 599	break;	×
600	}
601
602	Real alpha = 0.0;	193✔
603	Array xt(n), gt(n);	193✔
604	Real ft = 0.0;	193✔
605	if (!lineSearchWolfe(P, x, d, f, g, stpMax, alpha, xt, ft, gt)) {	193✔
NEW 606	ecType = EndCriteria::StationaryFunctionValue;	×
NEW 607	break;	×
608	}
609
610	// limited-memory update with the curvature safeguard
611	Array s = xt - x;	193✔
612	Array y = gt - g;	193✔
613	Real sy = DotProduct(s, y);	193✔
614	Real yy = DotProduct(y, y);	193✔
615
616	if (yy > 0.0 && sy > QL_EPSILON * yy) {	193✔
617	S.push_back(s);	193✔
618	Y.push_back(y);	193✔
619
620	if (S.size() > m_) {	193✔
621	S.pop_front();	136✔
622	Y.pop_front();	136✔
623	}
624	theta = yy / sy;	193✔
625	}
626
627	Real fOld = f;
628	x = xt;	193✔
629	f = ft;	193✔
630	g = gt;	193✔
631	++iter;	193✔
632
633	P.setCurrentValue(x);	386✔
634	P.setFunctionValue(f);
635
636	// relative function-reduction stop (SciPy's factr criterion)
637	Real denom = std::max(std::max(std::fabs(fOld), std::fabs(f)), Real(1.0));	193✔
638	if ((fOld - f) <= factr_ * QL_EPSILON * denom) {	193✔
639	ecType = EndCriteria::StationaryFunctionValue;	5✔
640	break;
641	}
642	}	193✔
643
644	P.setCurrentValue(x);	26✔
645	P.setFunctionValue(f);
646	P.setGradientNormValue(pgInf * pgInf);
647	return ecType;	13✔
648	}	13✔
649
650	}

lballabio / QuantLib / 28055112333

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous