Non-square Solve in PolynomialChaosExpansion

[asher-m]

This is another bug found by Eli. In particular, pce.PolynomialChoasExpansion.build calls into UncertainSCI.linalg.lstsq_loocv_error which, within this bug report, attempts to solve the problem $A x = b$ with $A$ of size $2 \times 35$ and $b$ of size $2 \times 2$ by calling numpy.linalg.solve. numpy.linalg.solve requires the last two dimensions of $A$ to be square (i.e., A.shape[-2] == A.shape[-1]) and throws an error.

MWE:

import numpy as np

import UncertainSCI.distributions as distributions
import UncertainSCI.pce as pce

samples = np.array(
    [[0.1, 0.4, 0.7],
     [0.2, 0.5, 0.8],
     [0.3, 0.6, 0.9]],
    dtype=np.float64
)
labels = ["param1", "param2", "param3"]

output = np.array(
    [[1.0, 4.0],
     [2.0, 5.0],
     [3.0, 6.0]],
    dtype=np.float64
)

o = 4
d = [
    distributions.BetaDistribution(alpha=1.0, beta=1.0, bounds=np.array((par.min(), par.max()), dtype=np.float64))
    for par in samples.T
]
p = pce.PolynomialChaosExpansion(distribution=d, order=o, plabels=labels)
p.set_samples(samples)
p.build(model_output=output)

This yields:

Console Log

File "C:\Users\asher\scm\uncertainsci-testfixes\UncertainSCI\pce.py", line 578, in build
  return self.build_pce_wlsq()
         ^^^^^^^^^^^^^^^^^^^^^
File "C:\Users\asher\scm\uncertainsci-testfixes\UncertainSCI\pce.py", line 309, in build_pce_wlsq
  self.accuracy_metrics['loocv'] = lstsq_loocv_error(V, self.model_output,
                                   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "C:\Users\asher\scm\uncertainsci-testfixes\UncertainSCI\utils\linalg.py", line 182, in lstsq_loocv_error
  A[m, P] @ np.linalg.solve(Rm, Qm.T @ bm))**2
            ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "C:\Users\asher\.miniconda3\envs\uncertainsci-testfixes\Lib\site-packages\numpy\linalg\_linalg.py", line 438, in solve
  _assert_stacked_square(a)
File "C:\Users\asher\.miniconda3\envs\uncertainsci-testfixes\Lib\site-packages\numpy\linalg\_linalg.py", line 246, in _assert_stacked_square
  raise LinAlgError('Last 2 dimensions of the array must be square')
numpy.linalg.LinAlgError: Last 2 dimensions of the array must be square

This may be as simple as doing some input checks on what is provided to the function by the user, or this may actually be a bug introduced by numpy2. I have not verified this with numpy < 2.