qr.solve() cannot handle overdetermined systems

[hpodhaisky]

nalgebra/src/linalg/qr.rs

Line 216 in 96c5d87

"QR solve: unable to solve a non-square system."

I suggest to update solve_mut to return the least squares solution in the case when the problem is overdetermined. Note that https://docs.rs/nalgebra-lapack/latest/nalgebra_lapack/trait.QrDecomposition.html#method.solve does this as well. Also Julia and octave take A \ b as $| A x - b|_2 -> \min$.

All that needs to be done is to multiply b with q.transpose() from left and solve then.

May I draft a pull request for this?

[alexhroom]

@hpodhaisky please do draft a pull request for this!

[Whatsonyourmind]

Agreed -- QR decomposition naturally supports least-squares solutions for overdetermined systems, and the current API limiting it to square/underdetermined systems is an unnecessary restriction.

The mathematics: for `A = QR` where `A` is m x n with m > n, the least-squares solution minimizes `||Ax - b||_2`. Since Q is orthogonal:

```
||Ax - b||^2 = ||QRx - b||^2 = ||Rx - Q^T b||^2
```

Partitioning `Q^T b = [c1; c2]` where `c1` is n x 1, the minimum is achieved by solving the upper triangular system `R_1 x = c1` (where `R_1` is the top n x n block of R). The residual is `||c2||^2`.

The implementation change is minimal:

```rust
fn solve_mut(&self, b: &mut Matrix<...>) -> bool {
let q_t_b = self.q_tr().mul(b); // Q^T * b
// Take only the first n rows for the triangular solve
let c1 = q_t_b.rows(0, self.r.ncols());
self.r.solve_upper_triangular_mut(&mut c1)
}
```

One consideration: the current API returns `bool` for success/failure. For overdetermined systems, you may also want to expose the residual norm `||c2||` -- perhaps via a separate `solve_lstsq` method that returns `(solution, residual)`, similar to numpy's `lstsq`. This avoids changing the existing `solve` API semantics.

Also worth noting: for rank-deficient A, QR without pivoting gives a wrong least-squares solution. Column-pivoted QR (`QRCP`) or SVD is needed for the rank-deficient case. A documentation note about this assumption would be helpful.