Fix psd solve vector rhs

[N-0-MAD]

Description

I fixed psd_solve_to_chol_solve so that it rewrites positive-definite Solve operations with a vector right-hand side (b_ndim=1) to use a Cholesky decomposition.

Previously, psd_solve_to_chol_solve only matched Solve operations with b_ndim=2. As a result, PyTensor left a vector solve such as:

A = pt.matrix("A", shape=(5, 5))
b = pt.vector("b", shape=(5,))
x = pt.linalg.solve(A, b, assume_a="pos", b_ndim=1)

as a Solve operation. This prevented the optimizer from reusing the Cholesky decomposition when the same matrix also appeared in a Cholesky-based expression such as a log-determinant.

I updated the rewrite to:

• Match both b_ndim=1 and b_ndim=2.
• Preserve the original b_ndim when constructing the triangular solves.
• Reuse the existing Cholesky/CholeskySolve machinery for vector RHS inputs without requiring users to reshape b to (n, 1).

I inspected the underlying SolveBase, CholeskySolve, and SolveTriangular implementations. They already support both b_ndim=1 and b_ndim=2.

I also tested the existing b[:, None] workaround. It successfully triggered the b_ndim=2 rewrite and produced a Cholesky-based graph.
I extended the existing test_psd_solve_with_chol test to parameterize over both b_ndim=1 and b_ndim=2, preserving the existing matrix coverage while adding regression coverage for vector RHS inputs.

Related Issue

• Closes psd_solve_to_chol_solve only tracks b_ndim=2, so solve(A, b, assume_a="pos") with a vector b factors A twice #2415
• Related to #

Checklist

• Checked that the pre-commit linting/style checks pass
• Included tests that prove the fix is effective or that the new feature works
• Added necessary documentation (docstrings and/or example notebooks)
• If you are a pro: each commit corresponds to a relevant logical change

Type of change

• Bug fix
• New feature / enhancement
• Documentation
• Maintenance
• Other (please specify):