[bodono]

Fixes #491.

Problem

proj_power_cone used absolute tolerances (POW_CONE_TOL = 1e-9) for the membership checks and the Newton stopping rule |f| < tol. With normalization and Anderson acceleration, the iterates that reach the projection in the issue's problem are tiny (~1e-6 to 1e-9). Newton stopped at r = 2.69e-9, f = -5.2e-10, which passes the absolute test but is ~24% from the root. That is exactly the reported x = 0.1238 vs. optimal 0.1. The residuals don't catch this because they assume the projection is exact, so SCS returned SCS_SOLVED, and it does so at every eps, not only 1e-8.

Fix (src/cones.c, proj_power_cone)

• Scale invariance: project v / ||v||_inf and rescale the result. Projection onto a cone is positively homogeneous, so the tolerances become relative to |v|.
• Cancellation-free x(r): when xh < 0, x(r) uses a form without cancellation, and the 1e-12 floor is gone. The floor overrode the true value when a component was tiny and negative.
• Bracketed Newton: phi(r) = x^a y^(1-a) is concave, so the sign of f = phi - r brackets the root. When a Newton step leaves [lo, hi] (from the left of a concave f it overshoots), an Illinois regula falsi step is used, then bisection. Bisection alone needed ~50 iterations when the root is within a few ulps of rh.
• Near the apex: right of the root, Newton on f(r)/r is tried (and used if it stays in the bracket). That form removes the trivial root at 0 and converges much faster than Newton on f, which is only linear near the double root. Restricting it to the near-apex regime saved a few projection evaluations on some inputs but made whole solves worse (see the solver sweep below), so it is taken whenever f < 0.
• Starting point: by default r = rh/2, with two better guesses where they apply.
  • Exactly one of xh, yh negative is the most common case for ADMM iterates. There, start at the root of the small-r approximation x(r) ~ a * rh * r / |xh|: r0 = (a * rh / |xh|)^(a / (1 - a)) * yh, or the mirror image. The true root is often tiny, and Newton from rh/2 took many shrinking steps to reach it.
  • For points just outside the cone (phi(rh) >= rh/2), such as converging iterates, start at r = rh. This only applies when |dx/dr| and |dy/dr| at rh are at most POW_CONE_START_SLOPE; otherwise the first step is too small to represent and looks converged.
  • Correctness doesn't depend on the guesses, since the bracket safeguards every step.
• Stopping rule: stop on the step, not on |f|. The step must be small relative to r and, scaled by sens = max(1, |dx/dr|, |dy/dr|), small relative to |v|. Near the apex x can be ~1000x more sensitive to r than r is. The relative part also keeps tiny components accurate relative to themselves when the root is near 0 (near the apex), which is the case that matters for the solves below. When the root is within a few ulps of rh and a component is tiny, no stopping rule in this r parametrization can do that; the error there is still small relative to |v| (the 2.3e-9 row below). The loop also stops when the bracket is that narrow ((hi - lo) * sens <= TOL) or only a few ulps wide.
• Constants: POW_CONE_TOL = 1e-12, POW_CONE_START_SLOPE = 1e3, POW_CONE_EPS = DBL_EPSILON (1e-6, 10, FLT_EPSILON under SFLOAT). POW_CONE_MAX_ITERS = 50.

Testing

• test_power_cone_small_scale uses the issue's exact data. It fails on master and passes here.

• test_power_cone_proj_tiny_component projects points just outside the cone that have one tiny positive component and p = 0.001 / 0.999, and checks them against 60-digit reference values. It fails on an intermediate commit of this PR that returned them unprojected.

• run_tests_direct and run_tests_indirect pass. The SFLOAT library builds with only a warning that is already on master.

• Issue reproducer: x = 0.09999999974, dist(s, K) = 1e-17 (was x = 0.1238, dist = 2.9e-3).

• Projection accuracy vs. a 60-digit decimal reference (KKT-verified on most sets; on the extreme-p / adversarial sets some true roots are ~1e-3000, so the reference is exact only to ~1e-59, far below double precision). Max error relative to ||v||_inf, mean Newton evaluations, and iteration-cap hits, 3000 points each:

  point set	master	this PR
  typical (`	v	` in 1e-3..1e3)
  just outside the cone	3.1e-5, 3.3 its, 77 at cap	1.2e-12, 1.7 its
  tiny positive component, just outside, p near 0/1	2.9e-2, 13.4 its, 732 at cap	2.3e-9, 4.5 its
  near the issue's dual solution (near-polar)	8.0e-8, 14.0 its	9.3e-13, 25.3 its
  extreme p (0.001, 0.01, 0.99, 0.999)	5.3e-5, 7.5 its, 363 at cap	2.6e-12, 5.4 its
  adversarial (badly scaled / near apex)	6.2e-1, 10.4 its, 479 at cap	1.9e-12, 10.7 its

  This PR never hits the 50-iteration cap in these sets.

• CVXPY benchmark: 52 models built with real 3-d power cones (approx=False), solved by scs-python built against master and against this PR on macOS arm64 / Accelerate, at eps = 1e-5 and 1e-8. The models are p-norm regression (p = 1.1, 1.5, 3), geometric-mean maximisation (uniform and skewed weights), power utility (p = 0.05, 0.3, 0.95), a power-penalty portfolio, CVXPY's pcp_1/2/3, and the issue's problem.

  • Accuracy: no objective is worse than master (vs. a tight Clarabel solve), and no model goes from solved to not solved. Two go the other way (skewed geometric means). The power-cone violation of the returned s is at most 0.06x the primal tolerance; on master it is up to 0.19x, and 2.4e4x on the issue's problem.
  • Speed (measured on 720f475, which differs from the final code only in when the f(r)/r step is taken): total solve time is dominated by iteration counts, which shift between builds in both directions because the ADMM trajectory is sensitive to last-bit differences in the projection. Per iteration, cone time is +8% vs. master (geometric mean over the 104 solves), and solve time is +8%. That ranges from 0.73-0.99x on p-norm (p = 1.1), portfolio and p = 0.95 utility models, to 1.17-1.45x on geometric-mean, pcp and p = 0.05/0.3 utility models. Before the last commit (speed-ups from profiling recorded projection inputs) cone time per iteration was +24%.
  • CVXPY's own SCS / power-cone tests (test_conic_solvers.py::TestSCS, test_pow_cone_nd.py, test_power_atom.py, test_power_tools.py, test_perspective.py, test_atoms.py) pass identically on both builds.

• Solver sweep, 390 problems: the 78 below with u jittered by 0-4%, eps = 1e-8. Not solved / total iterations: master (78-problem version) 61 / 6.2M; f(r)/r step only when Newton on f would more than halve r (720f475) 18 / 2.77M; only when both components are negative 26 / 3.48M; whenever f < 0 (final) 16 / 2.50M.

• Earlier solver sweep (measured on fd25f7a), min -z s.t. (u, 1, z) in K_p with u in 1e-12..1e-6 and p in {0.001, 0.01, 0.1, 0.9, 0.99, 0.999}, eps = 1e-8 (78 problems): master leaves 61 not solved (6.2M total iterations); this PR leaves 3 (476k). These problems are badly conditioned, so individual outcomes are sensitive to rounding-level changes in the projection. For p = 0.9 and small u, z* is below eps_abs, so errors in z there reflect the solver tolerance, not the projection.

• Reviewed five times by a second model; the last review approved the final code with no substantive findings. Earlier findings that are addressed: a converged step being discarded for bisection, cap hits when the root is near 0, a start guard that let tiny components through for small p, and the stopping-rule nits. One suggestion, stopping on the sensitivity-scaled step alone, made tiny-u solves worse, so it is combined with the relative test instead.

🤖 Generated with Claude Code