Non-identity DR Scaling
In this note we derive the update equations when using a non-identity scaling. Standard Douglas-Rachford splitting applied to the SCS problem is described in the Algorithm page. Now consider modifying DR splitting to use a diagonal matrix \(R\) instead of \(I\). This is useful because the matrix \(R\) can be selected to provide better convergence in practice. The algorithm then becomes:
which yields \(w^k \rightarrow u^\star + R^{-1} \mathcal{Q}(u^\star)\) where \(0 \in \mathcal{Q}(u^\star) + N_{\mathcal{C}_+}(u^\star)\).
This changes the first two steps of the procedure. The linear projection and the cone projection (explained next).
Cone projection
Note that \(R\) has to be chosen so that the cone projection can still be
computed. To do this we ensure that the entries of \(R > 0\) corresponding
to each cone are constant within that cone, with two exceptions. Both arise
only with the adaptive_diag_scale setting, which is
on by default and is described in Per-row refinement: each row of the
box cone can have its own entry, and the entries of a real PSD block can take
a rank-one form. (The zero and nonnegative cones are products of
one-dimensional cones, so each of their rows can have its own entry too.)
When the entries are constant within a cone, \(R\) does not affect the
projection onto that cone in any way, since for sub-cone \(\mathcal{K}\):
In other words, for these cones
For the two exceptions the projection does depend on \(R\), and it is computed exactly in the norm that \(R\) defines (see the Moreau decomposition below). The box cone projection solves a weighted problem in its scalar variable, and the PSD projection applies a diagonal congruence before and after the usual eigendecomposition.
Selecting \(R\)
For SCS we take
where \(I_n\) is \(n \times n\) identity, \(\rho_x \in \mathbf{R}\),
\(\rho_y \in \mathbf{R}^m\) and \(d \in \mathbf{R}\). The \(\rho_y\)
term includes the effect of the parameter scale, which is updated
heuristically to improve convergence. Basically
with two differences. First, each cone can modify this relationship slightly
(currently only the zero cone does, which sets
\(\rho_y = 1 / (1000\ \mathrm{scale})\)). Second, with
adaptive_diag_scale each row \(i\) has its own multiplier
\(f_i\), which acts like a scale for that row alone:
for a row outside the zero cone, and \(\rho_{y,i} = 1 / (1000\ \mathrm{scale}\ f_i)\) for a row of the zero cone. The multipliers start at one, so at the start of a solve on a new workspace they change nothing, and they are kept in \([10^{-3}, 10^{3}]\). How they are updated is described in Per-row refinement. So as scale increases, \(\rho_y\) decreases and vice-versa.
The quantity \(\rho_x\) is determined by the setting value
rho_x. The quantity \(\rho_y\) is determined by the setting value scale but is updated adaptively using the techniques
described in Dynamic scale updating. Finally, \(d\)
is determined by TAU_FACTOR in the code defined in glbopts.h.
Root-plus function
Finally, the root_plus function is modified to be the solution
of the following quadratic equation:
where \(R_{-1}\) corresponds to the first \(n+m\) entries of \(R\). Other than when computing \(\kappa\) (which does not affect the algorithm) this is the only place where \(d\) appears, so we have a lot of flexibility in how to choose it and it can even change from iteration to iteration. It is an open question on how best to select this parameter.
Moreau decomposition
The projection onto a cone under a diagonal scaling also satisfies a Moreau-style decomposition identity, as follows:
where \(\Pi_\mathcal{C}^R\) denotes the projection onto convex cone \(\mathcal{C}\) under the \(R\)-norm, which is defined as
This identity is useful when deriving cone projection routines, though most cones are either invariant to this or we enforce that \(R\) is constant across them. Note that the two components of the decomposition are \(R\)-orthogonal.
Dual vector
In order to get the dual vector \(v^k\) that contains \(s^k\) and \(\kappa^k\) we use:
and we have
by Moreau, and finally note that \(v^k \perp u^k\) from the fact that the Moreau decomposition is \(R\)-orthogonal.
Dynamic scale updating
The choice of the scale parameter can have a large impact on the
performance of the algorithm and the optimal choice is highly problem
dependent. SCS can dynamically adjust the scale parameter
on the fly via a heuristic procedure that can substantially improve convergence
in practice. This procedure is enabled by the adaptive_scale
setting. The procedure attempts to balance the convergence
rate of the primal residual with the dual residual. Loosely speaking, the
scale parameter will be increased if the primal residual is much larger
than the dual and decreased if the opposite is true.
The procedure does not run at every iteration. It runs at the iterations where
SCS computes the residuals to test for convergence, which is every 25
iterations starting at iteration 0 (CONVERGED_INTERVAL in
include/glbopts.h). At each of these iterations it takes one sample
of the imbalance between the two residuals and then decides whether to update.
Specifically, at such an iteration \(k\) consider the case where \(l\)
iterations have elapsed since the last update of the scale parameter
or of the per-row multipliers described below, and denote by
\((x, y, \tau) = u^k\) and \((0, s, \kappa) = v^k\).
Two measures of the imbalance between the primal and the dual residual are
blended. The first is how far each residual is from its own stopping
threshold, with \(\tau\) multiplied through:
so a value of one or less means that residual would pass the test. By default
we use the \(\ell_\infty\) norm for these quantities, which can be changed
using the SCALE_NORM constant in include/glbopts.h. The second
measure is the residuals themselves, \(\tilde r^k_p\) and
\(\tilde r^k_d\), in the space the iteration runs in (that is, after
equilibration when normalize is on), with
nothing divided out.
Balancing the distances \(d^k_p\) and \(d^k_d\) evens out how far each
side still has to go, which is the right target whenever both can be trusted.
A distance cannot be trusted when the residual is as large as the terms it is
measured against: such a residual has no correct digits, and its distance then
sits at a ceiling of about \(1 / \epsilon_\mathrm{rel}\) whatever the
iterate does. The clearest case is a feasibility problem, \(c = 0\) and
\(P = 0\), where the dual residual is \(A^\top y\) and is therefore
the only term it is measured against. Its distance is at the ceiling on every
iteration, so balancing the distances alone would read the dual as permanently
behind and drive the scale down to its lower limit. To handle this,
define the saturation of each side and the resulting weight
so that the scalar \(\theta^k\) is close to one when both residuals are small against their terms and zero when either is at its ceiling. The sample taken at iteration \(k\) is
Since the stopping thresholds enter the sample, the eps_abs and
eps_rel settings affect the updates of the scale, and
therefore the iterates, not only the point at which the solve stops. They
enter only through the ratio
\(\epsilon_\mathrm{abs} / \epsilon_\mathrm{rel}\): a common factor cancels
in \(d^k_p / d^k_d\) and in the saturations. So tightening both
tolerances by the same factor, for example from \(10^{-4}\) to
\(10^{-6}\) each, leaves the samples unchanged up to rounding, while
changing one of them alone changes the samples, and the solve can then follow
a different path from the first update on.
Now let \(S\) be the set of iterations at which a sample has been taken since the last update (or since the start of the solve), up to and including \(k\), and consider
ie, \(\beta\) is the geometric mean of the samples taken since the last
update. If this number is larger than 10 or smaller than 1/10 and if
sufficient iterations have passed since the last update (\(l\) at least
100, as determined by RESCALING_MIN_ITERS) then an update of the
scale parameter is triggered:
The presence of the square root is to prevent over-shooting the ‘optimal’
scale parameter, which could lead to oscillation. The limits are
MIN_SCALE_VALUE and MAX_SCALE_VALUE in
include/glbopts.h. If the scale is already at one of them and
\(\beta\) would push it further out, then no update of the scale
takes place (the per-row multipliers described below can still be updated on
their own). The limits apply to the updates only: the scale
setting can be any positive number.
With these values \(\beta\) is the geometric mean of four samples the first time an update is allowed after an earlier one (five for the first update of a solve, which also has the sample from iteration 0), and it gains one sample every 25 iterations for as long as no update is triggered.
Note that if the linear system is being solved using a
direct method, then every update requires a new factorization of the
perturbed matrix, so is somewhat expensive for larger problems and should be
done sparingly. Also, since changing the scale or the per-row
multipliers changes the operator we are using in DR splitting, every update
also resets the Anderson acceleration, whichever linear
solver is in use. After a reset the acceleration has to fill its memory again
before it produces a step: the first call after the reset only records the
iterate, and the first step comes acceleration_lookback calls later,
with one call every acceleration_interval iterations. With the default
settings that is 55 iterations without an accelerated step after each update,
and updates can be as close as 100 iterations apart.
The number of updates made in a solve is returned in scale_updates in
the info struct, and the final scale in scale.
Per-row refinement
The scale parameter above moves the whole \(R_y\) block together.
With adaptive_diag_scale each row \(i\) carries its own multiplier
\(f_i\) as well. It starts at one and divides that row’s entry of
\(\rho_y\), so row \(i\) behaves as if its scale were
\(f_i\ \mathrm{scale}\).
Both adaptive_scale and adaptive_diag_scale are on by default.
The multipliers are computed inside the procedure of the previous section,
which only runs when adaptive_scale is on, so
adaptive_diag_scale requires adaptive_scale. If
adaptive_scale is off then adaptive_diag_scale is switched off
as well, without a warning, and \(R\) stays fixed for the whole solve. If
adaptive_scale is on and adaptive_diag_scale is off then only
the scalar scale is updated.
The multipliers are updated from each row’s relative primal residual
by the damped multiplicative step
where \(g\) is the geometric mean of \(\hat r^k_{p,i}\) over all
\(m\) rows. A row whose relative residual is above the geometric mean
therefore gets a larger multiplier (a larger scale for that row), and a row
below it a smaller one. The result is clipped to \([10^{-3}, 10^{3}]\)
(the exponent and the range are DIAG_SCALE_DAMP,
DIAG_SCALE_MULT_MIN and DIAG_SCALE_MULT_MAX in
include/glbopts.h) and then constrained within each non-polyhedral
cone block to a form the cone projection can absorb. For most cone types that
is a single value, the arithmetic mean of the block’s clipped multipliers. For real PSD
blocks it is the richer rank-one family \(w_{ij} = \delta_i \delta_j\)
(a diagonal congruence \(X \to D X D\), which maps the cone to itself),
fitted per block by least squares in log space, with each \(\delta_i\)
clipped so that every \(w_{ij}\) stays in the same range. The PSD
projection under such a metric is computed exactly via a congruence sandwich
around the eigendecomposition. This is the dynamic analog of the row half of
Ruiz equilibration, using the residuals the algorithm actually produces rather
than the matrix norms.
The multipliers are updated whenever the scalar scale is. They are
also updated on their own, with a new factorization and a reset of the
acceleration exactly as for the scalar, when some multiplier would move by
more than the threshold used for the scalar (a factor of \(\sqrt{10}\))
and at least 100 iterations have passed since the last update of either
kind, so that a scalar scale sitting at its limit cannot freeze the
diagonal. The move is measured on the multipliers as they would be
installed, after each cone block has been given its form. A single row of a
second-order or PSD block that asks for a large step therefore does not cost a
factorization and a reset when the value its block ends up with has barely
changed. Rows of the zero, nonnegative and box cones each keep their own
multiplier, so for them nothing is averaged away. An update of the multipliers
on their own is counted in scale_updates (see Return information), but the
scale printed in the log and returned in the info struct does not
change.
The floor \(\epsilon_{\rm den}\) in the denominator is \(10^{-3}\)
(DEN_FLOOR_FRAC in include/glbopts.h) times the
root-mean-square over all rows of the denominator without the floor,
\(\max(|(Ax)_i|, |s_i|, |b_i \tau|)\). It is recomputed from the current
iterate each time, and the rows are taken in the space the iteration runs in
(that is, after equilibration when normalize is
on), like \(\tilde r^k_p\) above. It matters because a row
whose \(b_i\) is zero divides only by \(\max(|(Ax)_i|, |s_i|)\), and
both of those shrink as the iterates converge, so that row’s ratio inflates
even when it is converging perfectly well. Without the floor the profile
partly reports which terms a row happens to contain rather than how the row is
converging, and the metric chases that artifact. Rows whose denominators are
of ordinary size are unaffected.