1
\$\begingroup\$

I have implemented the Preconditioned Conjugate Gradient (PCG) method for solving a system of linear equations in Python and I would appreciate it if someone could verify its correctness since I am quite new to programming and using Python.

However, I am not entirely sure if it’s correct for all possible inputs. I would appreciate it if someone could review my code and provide feedback or point out any potential issues or have suggestions on how I can match it closer to the pseudocode image as it'll be easier for me to understand the code in the future by having it closely match the pseudocode.

import numpy as np
from scipy.sparse.linalg import spilu, LinearOperator
from scipy.sparse import csc_matrix
from scipy.linalg import solve_triangular

def PCG(A, b, L, U, tol, maxit):
    n = len(b)
    x = np.zeros(n)
    r = b - A @ x
    ILU1 = spilu(csc_matrix(A))
    M = LinearOperator((A.shape[0],A.shape[1]), ILU1.solve)
    z = M @ r
    d = z
    for i in range(maxit):
        Ad = A @ d
        alpha = np.dot(r.T, z) / np.dot(d.T, Ad)
        x = x + alpha * d
        r = r - alpha * Ad
        z_next = M @ r
        beta = np.dot(r.T, z_next) / np.dot(r.T, z)
        d = z_next + beta * d
        z = z_next
        if np.linalg.norm(r) < tol:
            break
    return x

Thank you in advance.

\$\endgroup\$
3
  • 1
    \$\begingroup\$ Code not implemented or not working as intended: Code Review is a community where programmers peer-review your working code to address issues such as security, maintainability, performance, and scalability. We require that the code be working correctly, to the best of the author's knowledge, before proceeding with a review. \$\endgroup\$ Commented Dec 9, 2023 at 13:11
  • 1
    \$\begingroup\$ If you've tested this and you think it works, please show your test results. \$\endgroup\$ Commented Dec 9, 2023 at 13:12
  • 2
    \$\begingroup\$ Can you provide a main function with examples so we can run it? It might seem a bit pedantic. But the thing is as a reviewer, I dont want to make assumptions on how you envision people using your code. Sometimes, the devil is in these details too. The way you think a piece of code will be used can be imoroved by reviewers here too. It also reduces the mental burden of the reviewer ensuring that you get a lot of replies which means a more comprehesive review. So, always post not just your code but a fully functional unit that can be copy pasted and compiled and run. \$\endgroup\$ Commented Dec 9, 2023 at 17:03

1 Answer 1

Reset to default
2
\$\begingroup\$

cite your reference

Thank you kindly for providing that .PNG algorithmic background as review context, it is very helpful. Besides the overall method, it explains why certain variable names were chosen.

The PCG function absolutely has to mention that original author, in the docstring. And it needs to spell out that we compute a preconditioned conjugate gradient. Which brings us to the next point...

condition number

Why are we conditioning at all? What is the domain that we are drawing our (possibly nearly) collinear matrices from, and how are we using them?

I would find a preconditioning argument much more convincing if we compute linalg.cond() both before and afterward, and verify the condition number was indeed reduced.

You didn't write a """docstring""", therefore the docstring isn't commenting at all on error sources nor on the consequences of producing "small" or "large" errors. As a user of this library, I would have no idea if my use case matches what this library addresses. And no idea if the library author had tested regions of the input- and output- space that I care about. Which brings us to the next issue...

test suite

You declined to provide a main() function that trivially exercises the target code, and we definitely don't see unit tests which exercise multiple regions of the input space. There is an opportunity to add small deltas to input matrices, and verify that your proposed method out-performs a baseline naïve method.

lint

Pep-8 asks that you spell it def pcg(..., lowercase. (And then there's the usual "math vs. snake_case" tension for a bunch of variables which we could more easily let slide, since local variables are not part of the Public API.)

optional type hinting

You know, given the empty docstring and the complex signature, we have arrived at the point where type annotation really isn't optional any more. Absent a citation in the code, a caller would have no idea what shape input parameters should be passed in, and what to expect as output. (One might reasonably expect the output matrix to be accompanied by a condition number or similar diagnostic, for example.)

I imagine that L, U were lower, upper. But I'm free to pass in anything I want; they are completely ignored and their meaning is obscure.

It is hard to tell if this codebase achieves its design goals, partly because a bunch of goals were left implicit rather than writing them down. You don't know if it properly works as intended, and neither do I. There are several techniques you could apply to improve our confidence in how it functions.

I would not be willing to delegate or accept maintenance tasks on this code.

\$\endgroup\$

You must log in to answer this question.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.