polysolve-python

Small Python binding for PolySolve's nonlinear and linear solver interfaces.

import numpy as np
import scipy.sparse
import polysolve


class Quadratic(polysolve.Problem):
    def value(self, x):
        y = x - np.array([-2.0, 3.0, 1.0])
        return float(y @ y)

    def gradient(self, x):
        return 2.0 * (x - np.array([-2.0, 3.0, 1.0]))

    def hessian(self, x):
        return 2.0 * scipy.sparse.eye(x.size, format="csc")


x, log = polysolve.minimize(
    Quadratic(),
    np.zeros(3),
    {
        "solver": "Newton",
        "line_search": {"method": "Backtracking"},
        "max_iterations": 100,
    },
    {"solver": "Eigen::SimplicialLDLT"},
)

print(x)
print(log)

Python subclasses must implement value(x), gradient(x), and hessian(x). Optional PolySolve callbacks such as solution_changed, stop, is_step_valid, and max_step_size can also be implemented on the subclass.

Linear solves

For a one-off linear system:

A = scipy.sparse.csc_matrix([[4.0, 1.0], [1.0, 3.0]])
b = np.array([1.0, 2.0])

x = polysolve.solve(A, b, {"solver": "Eigen::SimplicialLDLT"})

For repeated solves with the same matrix pattern:

solver = polysolve.LinearSolver({"solver": "Eigen::SimplicialLDLT"})
solver.analyze_pattern(A)
solver.factorise(A)  # factorize(A) is also available

x = solver.solve(b)
print(solver.info())