Add damp support to CG solver (closes #406) by gaoflow · Pull Request #752 · PyLops/pylops

gaoflow · 2026-05-29T10:19:53Z

Motivation

CG currently only solves $\mathbf{Op},\mathbf{x} = \mathbf{y}$. As discussed in #406, it should also be able to solve the damped system

$$(\mathbf{Op} + \epsilon\mathbf{I}),\mathbf{x} = \mathbf{y},$$

consistent with cgls/lsqr (which both expose damp). The damp argument used to exist on CG as an unused stub and was removed in v2; this PR re-introduces it with a working implementation.

This matches the issue's definition of done ("Added implementation with damp != 0") and the semantics @mrava87 confirmed there: solve $(A + \text{damp},I)x = y$, equivalent to adding damp*I to the operator and using damp=0.

Implementation

Damping is applied by adding $\epsilon\mathbf{c}$ to every application of the operator $\mathbf{Op},\mathbf{c}$ in the iterations (and $\epsilon\mathbf{x}_0$ to the initial residual when x0 is given) — i.e. the operator is implicitly replaced by $\mathbf{Op} + \epsilon\mathbf{I}$:

Opc = self.Op.matvec(self.c)
if self.damp != 0.0:
    Opc = Opc + self.damp * self.c

damp is threaded through cg(), CG.solve and CG.setup (default 0.0). With damp=0.0 the code path is identical to before, so existing behaviour is unchanged. Op must remain square and SPD, which $\mathbf{Op} + \epsilon\mathbf{I}$ is for $\epsilon \ge 0$.

Verification

Added test_cg_damp which checks, for both preallocate modes, the equivalence @mrava87 suggested in the issue:

internal damping cg(Op, y, damp=d) == adding d*I to the operator cg(Op + d*I, y, damp=0)
both == dense (Op + d*I)^{-1} y
damp=0 still recovers the undamped Op^{-1} y

Full pytests/test_solver.py passes (100 tests), ruff check clean.

The CG solver only solved Op x = y. Re-introduce the damp coefficient (removed in v2) so it can solve the damped system (Op + damp*I) x = y, consistent with cgls/lsqr. Damping is applied by adding damp*c to every operator application Op*c in the iterations (and to the initial residual when x0 is given); damp=0 reproduces the previous behaviour exactly. Thread damp through cg(), CG.solve and CG.setup, and document it. Add test_cg_damp verifying that internal damping matches both adding damp*I to the operator (with damp=0) and the dense (Op + damp*I)^-1 y, and that damp=0 still recovers the undamped solution.

codacy-production · 2026-05-29T10:20:50Z

Up to standards ✅

🟢 Issues 0 issues

Results:
0 new issues

View in Codacy

🟢 Metrics 0 complexity · 6 duplication

Metric Results

Complexity 0

Duplication 6

View in Codacy

_{NEW Get contextual insights on your PRs based on Codacy's metrics, along with PR and Jira context, without leaving GitHub. Enable AI reviewer}
_{TIP This summary will be updated as you push new changes.}

mrava87

@gaoflow thanks for a nice contribution :)

I left some minor comments, mostly to ensure consistency with other bits of codes in other solvers having the damp option.

mrava87 · 2026-05-29T20:55:49Z

    Solve the :math:`\mathbf{y} = \mathbf{Op}\,\mathbf{x}` problem using conjugate gradient
    iterations [1]_.

+    When a non-zero damping coefficient :math:`\epsilon` is provided to ``setup``/``solve``,


Follow the same style used in the notes of the other solvers. Simply:

Solve the :math:`\mathbf{y} = (\mathbf{Op} + + \epsilon\mathbf{I})\,\mathbf{x}` problem using conjugate gradient iterations [1]_, where :math:`\epsilon` is the damping coefficient.

Done in 6828c57 — reduced the CG Notes to the single damped-system sentence matching the other solvers.

mrava87 · 2026-05-29T20:56:58Z

        Initial guess
    niter : :obj:`int`, optional
        Number of iterations
+    damp : :obj:`float`, optional


Just Damping coefficient as in the other solvers

Done in 6828c57 — shortened to Damping coefficient.

mrava87 · 2026-05-29T20:57:21Z

        niter : :obj:`int`, optional
            Number of iterations (default to ``None`` in case a user wants to
            manually step over the solver)
+        damp : :obj:`float`, optional


Just Damping coefficient as in the other solvers

Done in 6828c57 — shortened to Damping coefficient.

mrava87 · 2026-05-29T21:01:14Z

                self.ncp.subtract(self.y, self.Op.matvec(x), out=self.r)
+            # account for the damping term in the initial residual
+            if self.damp != 0.0:
+                self.r -= self.damp * x


Follow the same pattern as self.r = self.y - self.Op.matvec(x) for not preallocate vs self.ncp.subtract(self.y, self.Op.matvec(x), out=self.r) for preallocate

Done in 6828c57 — setup now follows the non-preallocate self.r = self.r - damp * x vs preallocate self.ncp.subtract(..., out=self.r) pattern.

mrava87 · 2026-05-29T21:03:14Z


        """
        Opc = self.Op.matvec(self.c)
+        # solve the damped system (Op + damp*I) x = y by adding the damping


too long comment, just # add dampling contribution is enough

Done in 6828c57 — trimmed to # add damping contribution.

mrava87 · 2026-05-29T21:04:33Z

            internally as zero vector
        niter : :obj:`int`, optional
            Number of iterations
+        damp : :obj:`float`, optional


Same comment as above

Done in 6828c57 — trimmed to # add damping contribution.

mrava87 · 2026-05-29T21:05:34Z


+@pytest.mark.parametrize("par", [(par1), (par2)])
+def test_cg_damp(par):
+    """CG with damping solves the damped system (Op + damp*I) x = y (issue #406).


Remove (issue #406).

Done in 6828c57 — dropped the (issue #406). reference.

mrava87 · 2026-05-29T21:06:32Z

+    # adding damp*I to the operator and using damp=0 must match internal damping
+    Aop_damped = MatrixMult(A + damp * np.eye(n), dtype=par["dtype"])
+
+    for prealloc in [False, True]:


for preallocate in [False, True]: to be consistent with other tests naming

Done in 6828c57 — renamed the loop variable to for preallocate in [False, True]:.

mrava87 · 2026-05-29T21:07:51Z

+    A = np.random.normal(0, 1, (n, n)) + par["imag"] * np.random.normal(0, 1, (n, n))
+    A = np.conj(A).T @ A + np.eye(n)  # symmetric positive-definite
+    Aop = MatrixMult(A, dtype=par["dtype"])
+    y = np.random.normal(0, 1, n) + par["imag"] * np.random.normal(0, 1, n)


Like in other tests create x

x = np.ones(par["nx"]) + par["imag"] * np.ones(par["nx"]) y = Aop * x

Done in 6828c57 — now builds y = Aop * x from a known x.

- shorten damp docstrings to 'Damping coefficient' (cls_basic CG.setup/solve, basic.cg) and simplify CG Notes to a single damped-system sentence - follow preallocate vs non-preallocate residual-update pattern for the initial damping correction in setup - shorten the damping comment in step to '# add damping contribution' - test_cg_damp: build y from a known x (y = Aop * x), rename loop var to 'preallocate', drop issue ref from docstring

gaoflow · 2026-05-29T21:13:10Z

Thanks for the review! Addressed all points in 6828c57:

damp docstrings shortened to Damping coefficient in CG.setup, CG.solve and cg, and the CG Notes reduced to a single damped-system sentence.
Initial damping correction in setup now follows the non-preallocate (self.r = self.r - damp * x) vs preallocate (ncp.subtract(..., out=self.r)) pattern used elsewhere.
Trimmed the step comment to # add damping contribution.
test_cg_damp: now builds y = Aop * x from a known x, renamed the loop variable to preallocate, and dropped the issue ref from the docstring.

(The red build check is the self-hosted GPU runner failing on a ptxas version mismatch — Unsupported .version 8.8 — in unrelated test_kirchhoff/test_nonstatconvolve CUDA tests; all matrix, Lint and Azure builds are green.)

mrava87 · 2026-05-30T07:35:27Z

Thanks for the review! Addressed all points in 6828c57:

damp docstrings shortened to Damping coefficient in CG.setup, CG.solve and cg, and the CG Notes reduced to a single damped-system sentence.

Initial damping correction in setup now follows the non-preallocate (self.r = self.r - damp * x) vs preallocate (ncp.subtract(..., out=self.r)) pattern used elsewhere.

Trimmed the step comment to # add damping contribution.

test_cg_damp: now builds y = Aop * x from a known x, renamed the loop variable to preallocate, and dropped the issue ref from the docstring.

(The red build check is the self-hosted GPU runner failing on a ptxas version mismatch — Unsupported .version 8.8 — in unrelated test_kirchhoff/test_nonstatconvolve CUDA tests; all matrix, Lint and Azure builds are green.)

Thanks. Please reply to each comment with link to commit you fixed the issue.

mrava87 requested changes May 29, 2026

View reviewed changes

Conversation

gaoflow commented May 29, 2026

Motivation

Implementation

Verification

Uh oh!

codacy-production Bot commented May 29, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Up to standards ✅

Uh oh!

mrava87 left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

gaoflow commented May 29, 2026

Uh oh!

mrava87 commented May 30, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

codacy-production Bot commented May 29, 2026 •

edited

Loading