Geometry-aware point removal bug

[lindonroberts]

After adding xnew, we have p+2 points in the model. Then selecting the best point to remove based on poisedness ends up calculating linear Lagrange polynomials in a p+1 dimensional subspace rather than the proper p dimensional subspace. This happens because of numerical errors, so xnew doesn't appear to be exactly in the existing space, and model.R ends up extremely ill-conditioned (if singular we revert to distance-only point removal).

Need to update this process, perhaps using regression Lagrange polynomials in the original p dimensional subspace for the p+2 points? Important to check this doesn't decrease overall performance of the code.

Code snippet making this clear:

import numpy as np
from dfbgn.model import InterpSet
from dfbgn.hessian import Hessian
from dfbgn.trust_region import trsbox

def objfun(x):  # f(x) = ||x||^2, simple least-squares
    return x.copy()

np.random.seed(0)
n = 10
fixed_block = 3
delta = 0.2
maxfun = 100
args = ()
scaling_changes = None
x0 = np.ones((n,), dtype=float)

r0 = objfun(x0)
nf = 1
xl = -1e20 * np.ones((n,))  # unconstrained
xu = 1e20 * np.ones((n,))  # unconstrained

# Start of dfbgn.solve_main(...)
model = InterpSet(fixed_block, x0, r0, xl, xu)
exit_info, nf = model.initialise_interp_set(delta, objfun, args, scaling_changes, nf, maxfun)
print("Initially, model has %g directions" % model.p)

xk = model.xopt()
fk = model.fopt()
interp_ok, ck, Jk = model.interpolate_mini_models(jac_in_full_space=False)
gk = 2.0 * Jk.T.dot(ck)
Hk = Hessian(model.p, vals=2.0 * np.dot(Jk.T, Jk))
sk_red, _, crvmin = trsbox(np.zeros((model.p,)), gk, Hk, -1e20 * np.ones((model.p,)), 1e20 * np.ones((model.p,)), delta)
sk_full = model.project_to_full_space(sk_red)

xnew = xk + sk_full
nf += 1
rnew = objfun(xnew)
fnew = np.dot(rnew, rnew)

pred_reduction = -np.dot(sk_red, gk + 0.5 * Hk.vec_mul(sk_red))
actual_reduction = fk - fnew
ratio = actual_reduction / pred_reduction

print("fk = %g, fnew = %g, ratio = %g" % (fk, fnew, ratio))

model.append_point(xnew, rnew)  # updates xopt
xnew_appended = True

min_npt_to_drop = (1 if xnew_appended else 0) + (2 if fixed_block is None else 1)
alpha = 1.0  # on successful steps
ndirs_to_keep = min(int(alpha * model.p), model.p - min_npt_to_drop)
ndirs_to_keep = max(0, ndirs_to_keep)
ndirs_to_drop = model.p - ndirs_to_keep

print("After adding xnew, model has %g directions, dropping %g" % (model.p, ndirs_to_drop))

try:
    sigmas = model.poisedness_of_each_point(delta=delta)
    print("sigmas =", sigmas)
    print("abs(Rdiag) =", np.abs(np.diag(model.R)))
    sqdists = np.square(model.distances_to_xopt(include_kopt=True))  # ||yt-xk||^2
    vals = sigmas * np.maximum(sqdists**2 / delta**4, 1)  # BOBYQA point to remove criterion
    vals[model.kopt] = -1.0  # make sure kopt is never selected
    print("vals =", vals)
except np.linalg.LinAlgError:
    print("Poisedness calculation failed, using distance to xopt")
    # If poisedness calculation fails, revert to dropping furthest points
    vals = np.square(model.distances_to_xopt(include_kopt=True))  # ||yt-xk||^2
    vals[model.kopt] = -1.0  # make sure kopt is never selected
    print("vals =", vals)
for i in range(ndirs_to_drop):
    k = np.argmax(vals)
    vals = np.delete(vals, k)  # keep vals indices in line with indices of model.points
    model.remove_point(k, check_not_kopt=True)

print("Finish iteration with %g directions" % model.p)
print("At start of next iteration, will add new directions until have %g directions" % fixed_block)