Decompose Tridiagonal Solve into core steps

Author: ricardoV94

pytensor/tensor/_linalg/solve/rewriting.py

@@ -4,6 +4,10 @@
 from pytensor.graph.rewriting.basic import copy_stack_trace, in2out, node_rewriter
 from pytensor.scan.op import Scan
 from pytensor.scan.rewriting import scan_seqopt1
+from pytensor.tensor._linalg.solve.tridiagonal import (
+    tridiagonal_lu_factor,
+    tridiagonal_lu_solve,
+)
 from pytensor.tensor.basic import atleast_Nd
 from pytensor.tensor.blockwise import Blockwise
 from pytensor.tensor.elemwise import DimShuffle
@@ -16,13 +20,17 @@
 def decompose_A(A, assume_a):
     if assume_a == "gen":
         return lu_factor(A, check_finite=False)
+    elif assume_a == "tridiagonal":
+        return tridiagonal_lu_factor(A)
     else:
         raise NotImplementedError
 
 
 def solve_lu_decomposed_system(A_decomp, b, b_ndim, assume_a, transposed=False):
     if assume_a == "gen":
         return lu_solve(A_decomp, b, b_ndim=b_ndim, trans=transposed)
+    elif assume_a == "tridiagonal":
+        return tridiagonal_lu_solve(A_decomp, b, b_ndim=b_ndim, transposed=transposed)
     else:
         raise NotImplementedError
 
@@ -58,7 +66,7 @@ def find_solve_clients(var, assume_a):
 
     assume_a = node.op.core_op.assume_a
 
-    if assume_a != "gen":
+    if assume_a not in {"gen", "tridiagonal"}:
         return None
 
     A, _ = get_root_A(node.inputs[0])
@@ -139,7 +147,7 @@ def scan_pushout_solve_lu_decomposition(fgraph, node):
             if (
                 isinstance(inner_node.op, Blockwise)
                 and isinstance(inner_node.op.core_op, Solve)
-                and inner_node.op.core_op.assume_a == "gen"
+                and inner_node.op.core_op.assume_a not in {"gen", "tridiagonal"}
             ):
                 # TODO: Move transpose from graph to Low level solve Op
                 A, _ = inner_node.inputs

pytensor/tensor/_linalg/solve/tridiagonal.py

@@ -0,0 +1,154 @@
+import numpy as np
+from scipy.linalg import get_lapack_funcs
+
+from pytensor.graph import Apply, Op
+from pytensor.tensor.basic import as_tensor, diagonal
+from pytensor.tensor.blockwise import Blockwise
+from pytensor.tensor.type import tensor, vector
+
+
+class LUFactorTridiagonal(Op):
+    """Compute LU factorization of a tridiagonal matrix (lapack gttrf)"""
+
+    __props__ = (
+        "overwrite_dl",
+        "overwrite_d",
+        "overwrite_du",
+    )
+    gufunc_signature = "(dl),(d),(dl)->(dl),(d),(dl),(du2),(d)"
+
+    def __init__(self, overwrite_dl=False, overwrite_d=False, overwrite_du=False):
+        self.destroy_map = dm = {}
+        if overwrite_dl:
+            dm[0] = [0]
+        if overwrite_d:
+            dm[1] = [1]
+        if overwrite_du:
+            dm[2] = [2]
+        self.overwrite_dl = overwrite_dl
+        self.overwrite_d = overwrite_d
+        self.overwrite_du = overwrite_du
+        super().__init__()
+
+    def make_node(self, dl, d, du):
+        dl, d, du = map(as_tensor, (dl, d, du))
+
+        if not all(inp.type.ndim == 1 for inp in (dl, d, du)):
+            raise ValueError("Diagonals must be vectors")
+
+        ndl, nd, ndu = (inp.type.shape[-1] for inp in (dl, d, du))
+        n = (
+            ndl + 1
+            if ndl is not None
+            else (nd if nd is not None else (ndu + 1 if ndu is not None else None))
+        )
+        dummy_arrays = [np.zeros((), dtype=inp.type.dtype) for inp in (dl, d, du)]
+        out_dtype = get_lapack_funcs("gttrf", dummy_arrays).dtype
+        outputs = [
+            vector(shape=(None if n is None else (n - 1),), dtype=out_dtype),
+            vector(shape=(n,), dtype=out_dtype),
+            vector(shape=(None if n is None else n - 1,), dtype=out_dtype),
+            vector(shape=(None if n is None else n - 2,), dtype=out_dtype),
+            vector(shape=(n,), dtype=np.int32),
+        ]
+        return Apply(self, [dl, d, du], outputs)
+
+    def perform(self, node, inputs, output_storage):
+        gttrf = get_lapack_funcs("gttrf", dtype=node.outputs[0].type.dtype)
+        dl, d, du, du2, ipiv, _ = gttrf(
+            *inputs,
+            overwrite_dl=self.overwrite_dl,
+            overwrite_d=self.overwrite_d,
+            overwrite_du=self.overwrite_du,
+        )
+        output_storage[0][0] = dl
+        output_storage[1][0] = d
+        output_storage[2][0] = du
+        output_storage[3][0] = du2
+        output_storage[4][0] = ipiv
+
+
+class SolveLUFactorTridiagonal(Op):
+    """Solve a system of linear equations with a tridiagonal coefficient matrix (lapack gttrs)."""
+
+    __props__ = ("b_ndim", "overwrite_b", "transposed")
+
+    def __init__(self, b_ndim: int, transposed: bool, overwrite_b=False):
+        if b_ndim not in (1, 2):
+            raise ValueError("b_ndim must be 1 or 2")
+        if b_ndim == 1:
+            self.gufunc_signature = "(dl),(d),(dl),(du2),(d),(d)->(d)"
+        else:
+            self.gufunc_signature = "(dl),(d),(dl),(du2),(d),(d,rhs)->(d,rhs)"
+        if overwrite_b:
+            self.destroy_map = {0: [5]}
+        self.b_ndim = b_ndim
+        self.transposed = transposed
+        self.overwrite_b = overwrite_b
+        super().__init__()
+
+    def make_node(self, dl, d, du, du2, ipiv, b):
+        dl, d, du, du2, ipiv, b = map(as_tensor, (dl, d, du, du2, ipiv, b))
+
+        if b.type.ndim != self.b_ndim:
+            raise ValueError("Wrang number of dimensions for input b.")
+
+        if not all(inp.type.ndim == 1 for inp in (dl, d, du, du2, ipiv)):
+            raise ValueError("Inputs must be vectors")
+
+        ndl, nd, ndu, ndu2, nipiv = (
+            inp.type.shape[-1] for inp in (dl, d, du, du2, ipiv)
+        )
+        nb = b.type.shape[0]
+        n = (
+            ndl + 1
+            if ndl is not None
+            else (
+                nd
+                if nd is not None
+                else (
+                    ndu + 1
+                    if ndu is not None
+                    else (
+                        ndu2 + 2
+                        if ndu2 is not None
+                        else (nipiv if nipiv is not None else nb)
+                    )
+                )
+            )
+        )
+        dummy_arrays = [
+            np.zeros((), dtype=inp.type.dtype) for inp in (dl, d, du, du2, ipiv)
+        ]
+        # Seems to always be float64?
+        out_dtype = get_lapack_funcs("gttrs", dummy_arrays).dtype
+        if self.b_ndim == 1:
+            output_shape = (n,)
+        else:
+            output_shape = (n, b.type.shape[-1])
+
+        outputs = [tensor(shape=output_shape, dtype=out_dtype)]
+        return Apply(self, [dl, d, du, du2, ipiv, b], outputs)
+
+    def perform(self, node, inputs, output_storage):
+        gttrs = get_lapack_funcs("gttrs", dtype=node.outputs[0].type.dtype)
+        x, _ = gttrs(
+            *inputs,
+            overwrite_b=self.overwrite_b,
+            trans="N" if not self.transposed else "T",
+        )
+        output_storage[0][0] = x
+
+
+def tridiagonal_lu_factor(a):
+    # Return the decomposition of A implied by a solve tridiagonal
+    dl, d, du = (diagonal(a, offset=o, axis1=-2, axis2=-1) for o in (-1, 0, 1))
+    dl, d, du, du2, ipiv = Blockwise(LUFactorTridiagonal())(dl, d, du)
+    return dl, d, du, du2, ipiv
+
+
+def tridiagonal_lu_solve(a_diagonals, b, *, b_ndim: int, transposed: bool = False):
+    dl, d, du, du2, ipiv = a_diagonals
+    return Blockwise(SolveLUFactorTridiagonal(b_ndim=b_ndim, transposed=transposed))(
+        dl, d, du, du2, ipiv, b
+    )

pytensor/tensor/rewriting/subtensor.py

@@ -20,6 +20,7 @@
 from pytensor.scalar import constant as scalar_constant
 from pytensor.tensor.basic import (
     Alloc,
+    ExtractDiag,
     Join,
     MakeVector,
     ScalarFromTensor,
@@ -28,6 +29,7 @@
     as_tensor,
     cast,
     concatenate,
+    full,
     get_scalar_constant_value,
     get_underlying_scalar_constant_value,
     register_infer_shape,
@@ -2163,3 +2165,82 @@ def ravel_multidimensional_int_idx(fgraph, node):
     "numba",
     use_db_name_as_tag=False,  # Not included if only "specialize" is requested
 )
+
+
+@register_canonicalize
+@register_stabilize
+@register_specialize
+@node_rewriter([ExtractDiag])
+def extract_diag_of_diagonal_set_subtensor(fgraph, node):
+    def is_contant_arange(var) -> bool:
+        if not (isinstance(var, TensorConstant) and var.type.ndim == 1):
+            return False
+
+        data = var.data
+        start, stop = data[0], data[-1] + 1
+        return data.size == (stop - start) and (data == np.arange(start, stop)).all()
+
+    [diag_x] = node.inputs
+    if not (
+        diag_x.owner is not None
+        and isinstance(diag_x.owner.op, AdvancedIncSubtensor)
+        and diag_x.owner.op.set_instead_of_inc
+    ):
+        return None
+
+    x, y, *idxs = diag_x.owner.inputs
+
+    if not (
+        x.type.ndim >= 2
+        and None not in x.type.shape[-2:]
+        and x.type.shape[-2] == x.type.shape[-1]
+    ):
+        # For now we only support rewrite with static square shape for x
+        return None
+
+    op = node.op
+    if op.axis2 > len(idxs):
+        return None
+
+    # Check all non-axis indices are full slices
+    axis = {op.axis1, op.axis2}
+    if not all(is_full_slice(idx) for i, idx in enumerate(idxs) if i not in axis):
+        return None
+
+    # Check axis indices are arange we would expect from setting on the diagonal
+    axis1_idx = idxs[op.axis1]
+    axis2_idx = idxs[op.axis2]
+    if not (is_contant_arange(axis1_idx) and is_contant_arange(axis2_idx)):
+        return None
+
+    dim_length = x.type.shape[-1]
+    offset = op.offset
+    start_stop1 = (axis1_idx.data[0], axis1_idx.data[-1] + 1)
+    start_stop2 = (axis2_idx.data[0], axis2_idx.data[-1] + 1)
+    orig_start1, orig_start2 = start_stop1[0], start_stop2[0]
+
+    if offset < 0:
+        # The logic for checking if we are selecting or not a diagonal for negative offset is the same
+        # as the one with positive offset but swapped axis
+        start_stop1, start_stop2 = start_stop2, start_stop1
+        offset = -offset
+
+    start1, stop1 = start_stop1
+    start2, stop2 = start_stop2
+    if (
+        start1 == 0
+        and start2 == offset
+        and stop1 == dim_length - offset
+        and stop2 == dim_length
+    ):
+        # We are extracting the just written diagonal
+        if y.type.ndim == 0 or y.type.shape[-1] == 1:
+            # We may need to broadcast y
+            y = full((*x.shape[:-2], dim_length - offset), y, dtype=x.type.dtype)
+        return [y]
+    elif (orig_start2 - orig_start1) != op.offset:
+        # Some other diagonal was written, ignore it
+        return [op(x)]
+    else:
+        # A portion, but no the whole diagonal was written, don't do anything
+        return None

pytensor/tensor/subtensor.py

@@ -3021,12 +3021,7 @@ def make_node(self, x, y, *inputs):
         return Apply(
             self,
             (x, y, *new_inputs),
-            [
-                tensor(
-                    dtype=x.type.dtype,
-                    shape=tuple(1 if s == 1 else None for s in x.type.shape),
-                )
-            ],
+            [x.type()],
         )
 
     def perform(self, node, inputs, out_):

tests/tensor/rewriting/test_subtensor.py

@@ -1,3 +1,5 @@
+import random
+
 import numpy as np
 import pytest
 
@@ -2415,3 +2417,37 @@ def test_unknown_step(self):
             f(test_x, -2),
             test_x[0:3:-2, -1:-6:2, ::],
         )
+
+
+def test_extract_diag_of_diagonal_set_subtensor():
+    A = pt.full((2, 6, 6), np.nan)
+    rows = pt.arange(A.shape[-2])
+    cols = pt.arange(A.shape[-1])
+    write_offsets = [-2, -1, 0, 1, 2]
+    # Randomize order of write operations, to make sure rewrite is not sensitive to it
+    random.shuffle(write_offsets)
+    for offset in write_offsets:
+        value = offset + 0.1 * offset
+        if offset == 0:
+            A = A[..., rows, cols].set(value)
+        elif offset > 0:
+            A = A[..., rows[:-offset], cols[offset:]].set(value)
+        else:
+            offset = -offset
+            A = A[..., rows[offset:], cols[:-offset]].set(value)
+    # Add a partial diagonal along offset 3
+    A = A[..., rows[1:-3], cols[4:]].set(np.pi)
+
+    read_offsets = [-2, -1, 0, 1, 2, 3]
+    outs = [A.diagonal(offset=offset, axis1=-2, axis2=-1) for offset in read_offsets]
+    rewritten_outs = rewrite_graph(outs, include=("ShapeOpt", "canonicalize"))
+
+    # Every output should just be an Alloc with value
+    expected_outs = []
+    for offset in read_offsets[:-1]:
+        value = np.asarray(offset + 0.1 * offset, dtype=A.type.dtype)
+        expected_outs.append(pt.full((np.int64(2), np.int8(6 - abs(offset))), value))
+    # The partial diagonal shouldn't be rewritten
+    expected_outs.append(outs[-1])
+
+    assert equal_computations(rewritten_outs, expected_outs)