Repository: gonum/gonum (language: go)

Issue:
mat: calling qr.Factorize leads to OOM for matrixes with many rows
### What are you trying to do?

After upgrading gonum to 0.15.0 calls to `(VecDense).SolveVec(a,b)` where `a.Dims() == (10k+, 3)`  leads to OOM.
The stack where OOM happens:
```
runtime.systemstack_switch()
        /nix/store/8yw3g52r95h7cv09lcrran92n212997b-go-1.22.0/share/go/src/runtime/asm_amd64.s:474 +0x8 fp=0xc000087230 sp=0xc000087220 pc=0x4788a8
runtime.(*mheap).alloc(0x3667b52000?, 0x1b33da9?, 0x80?)
        /nix/store/8yw3g52r95h7cv09lcrran92n212997b-go-1.22.0/share/go/src/runtime/mheap.go:958 +0x5b fp=0xc000087278 sp=0xc000087230 pc=0x42fc1b
runtime.(*mcache).allocLarge(0x100?, 0x3667b50b88, 0x1?)
        /nix/store/8yw3g52r95h7cv09lcrran92n212997b-go-1.22.0/share/go/src/runtime/mcache.go:234 +0x85 fp=0xc0000872c0 sp=0xc000087278 pc=0x41ce45
runtime.mallocgc(0x3667b50b88, 0xd67320, 0x1)
        /nix/store/8yw3g52r95h7cv09lcrran92n212997b-go-1.22.0/share/go/src/runtime/malloc.go:1165 +0x597 fp=0xc000087348 sp=0xc0000872c0 pc=0x4141d7
runtime.makeslice(0xc0000124b0?, 0x18?, 0xc0000874b0?)
        /nix/store/8yw3g52r95h7cv09lcrran92n212997b-go-1.22.0/share/go/src/runtime/slice.go:107 +0x49 fp=0xc000087370 sp=0xc000087348 pc=0x459989
gonum.org/v1/gonum/mat.NewDense(...)
        /home/ericwenn/go/pkg/mod/gonum.org/v1/gonum@v0.15.1-0.20240413203616-1b7d9ca04ac9/mat/dense.go:60
gonum.org/v1/gonum/mat.(*QR).updateQ(0xc000087630)
        /home/ericwenn/go/pkg/mod/gonum.org/v1/gonum@v0.15.1-0.20240413203616-1b7d9ca04ac9/mat/qr.go:107 +0x78 fp=0xc0000874c0 sp=0xc000087370 pc=0xb5fd38
gonum.org/v1/gonum/mat.(*QR).factorize(0xc000087630, {0x1070598, 0xc0004c6540}, 0x49)
        /home/ericwenn/go/pkg/mod/gonum.org/v1/gonum@v0.15.1-0.20240413203616-1b7d9ca04ac9/mat/qr.go:101 +0x316 fp=0xc0000875e0 sp=0xc0000874c0 pc=0xb5fc36
gonum.org/v1/gonum/mat.(*QR).Factorize(...)
        /home/ericwenn/go/pkg/mod/gonum.org/v1/gonum@v0.15.1-0.20240413203616-1b7d9ca04ac9/mat/qr.go:82
gonum.org/v1/gonum/mat.(*Dense).Solve(0xc0004c6700, {0x1070598, 0xc0004c6540}, {0x1070598, 0xc0004c6740})
        /home/ericwenn/go/pkg/mod/gonum.org/v1/gonum@v0.15.1-0.20240413203616-1b7d9ca04ac9/mat/solve.go:65 +0x245 fp=0xc000087700 sp=0xc0000875e0 pc=0xb63545
gonum.org/v1/gonum/mat.(*VecDense).SolveVec(0xc000110360, {0x1070598, 0xc0004c6540}, {0x10735c0, 0xc0001103c0})
        /home/ericwenn/go/pkg/mod/gonum.org/v1/gonum@v0.15.1-0.20240413203616-1b7d9ca04ac9/mat/solve.go:118 +0x6fa fp=0xc000087930 sp=0xc000087700 pc=0xb63e9a
``` 
I've traced this down to changes in `mat.QR.updateQ` [here](https://github.com/gonum/gonum/commit/45b74210d616dd61bb893f6e8e4af13879b1d6bb#diff-5bf82974a09725d659ece49d56eff0b81eafc1a476a618949ec743544e329aaeR104-R117).

What `updateQ` does is allocate a new `Dense` with dimensions `MxM`, which for our value on M would mean allocating 100M cell matrix. Based on my understanding of the code, the resulting matrix should have dimensions (`MxN`).


### What version of Go and Gonum are you using?
<!--
Paste the output of `go version` and if you are installing Gonum from source, paste
the output of `(cd $(go env GOPATH)/src/gonum.org/v1/gonum && git rev-parse HEAD)`.
If you are using modules, also paste the output of `grep gonum.org/v1/gonum go.sum`,
executed in the root of your dependent module.
-->

```
$ go version
go version go1.22.0 linux/amd64

$ grep gonum.org/v1/gonum go.sum
gonum.org/v1/gonum v0.15.0 h1:2lYxjRbTYyxkJxlhC+LvJIx3SsANPdRybu1tGj9/OrQ=
gonum.org/v1/gonum v0.15.0/go.mod h1:xzZVBJBtS+Mz4q0Yl2LJTk+OxOg4jiXZ7qBoM0uISGo=
```

### Does this issue reproduce with the current master?

Yes



Discussion / hints:
Can you provide a minimal reproducer for testing and fix the link in the report please?

Note also from the documentation that Q is m×m.
I believe the issue arises when one tries to solve very large systems.

The issue is introduced by [this commit](https://github.com/gonum/gonum/commit/45b74210d616dd61bb893f6e8e4af13879b1d6bb) which explicitly computes and stores the `Q` matrix (instead of just the reflectors `tau`).

The `Q` matrix doesn't seem to be explicitly needed in the QR `SolveVec`, but the new update forces its computation, which in practice can limit the size of the linear systems the function can solve.
Minimal reproduction test:
```go
func TestReproduction(t *testing.T) {
	const (
		M = 100_000
		N = 3
	)
	systemMat := mat.NewDense(M, N, make([]float64, M*N))
	measurementsVec := mat.NewVecDense(M, make([]float64, M))

	var results mat.VecDense
	err := results.SolveVec(systemMat, measurementsVec)

	fmt.Println(err)
}
```

On version `0.14.0` and below this correctly returns error `matrix singular or near-singular with condition number +Inf`, on version `0.15.0` it crashes due to OOM.

PS I've updated the link in original issue, thanks for pointing that out.
Thanks. Yeah, for situations where the Q is not needed, this is a problem. I think we can get the performance benefit of the single pre-calculation by lazily calculating, but storing the value so the work is only done once. I'll take a look in the weekend.
> Thanks. Yeah, for situations where the Q is not needed, this is a problem. I think we can get the performance benefit of the single pre-calculation by lazily calculating, but storing the value so the work is only done once. I'll take a look in the weekend.

That would be great, thanks for the quick response!

Fix patch:
diff --git a/AUTHORS b/AUTHORS
index 791b69856..8e5896dbf 100644
--- a/AUTHORS
+++ b/AUTHORS
@@ -125,6 +125,7 @@ The University of Washington
 Thomas Berg <tomfuture@gmail.com>
 Tobin Harding <me@tobin.cc>
 Tom Payne <twpayne@gmail.com>
+Tristan Nicholls <tvk.nicholls@gmail.com>
 Valentin Deleplace <deleplace2015@gmail.com>
 Vincent Thiery <vjmthiery@gmail.com>
 Vladimír Chalupecký <vladimir.chalupecky@gmail.com>
diff --git a/CONTRIBUTORS b/CONTRIBUTORS
index a7e62e1e4..e367595b9 100644
--- a/CONTRIBUTORS
+++ b/CONTRIBUTORS
@@ -128,6 +128,7 @@ Tamir Hyman <hyman.tamir@gmail.com>
 Thomas Berg <tomfuture@gmail.com>
 Tobin Harding <me@tobin.cc>
 Tom Payne <twpayne@gmail.com>
+Tristan Nicholls <tvk.nicholls@gmail.com>
 Valentin Deleplace <deleplace2015@gmail.com>
 Vincent Thiery <vjmthiery@gmail.com>
 Vladimír Chalupecký <vladimir.chalupecky@gmail.com>
diff --git a/mat/qr.go b/mat/qr.go
index af99dbcaa..7f8fec8f6 100644
--- a/mat/qr.go
+++ b/mat/qr.go
@@ -31,8 +31,13 @@ func (qr *QR) Dims() (r, c int) {
 	return qr.qr.Dims()
 }
 
-// At returns the element at row i, column j.
+// At returns the element at row i, column j. At will panic if the receiver
+// does not contain a successful factorization.
 func (qr *QR) At(i, j int) float64 {
+	if !qr.isValid() {
+		panic(badQR)
+	}
+
 	m, n := qr.Dims()
 	if uint(i) >= uint(m) {
 		panic(ErrRowAccess)
@@ -41,6 +46,20 @@ func (qr *QR) At(i, j int) float64 {
 		panic(ErrColAccess)
 	}
 
+	if qr.q == nil || qr.q.IsEmpty() {
+		// Calculate Qi, Q i-th row
+		qi := getFloat64s(m, true)
+		qr.qRowTo(i, qi)
+
+		// Compute QR(i,j)
+		var val float64
+		for k := 0; k <= j; k++ {
+			val += qi[k] * qr.qr.at(k, j)
+		}
+		putFloat64s(qi)
+		return val
+	}
+
 	var val float64
 	for k := 0; k <= j; k++ {
 		val += qr.q.at(i, k) * qr.qr.at(k, j)
@@ -48,6 +67,25 @@ func (qr *QR) At(i, j int) float64 {
 	return val
 }
 
+// qRowTo extracts the i-th row of the orthonormal matrix Q from a QR
+// decomposition.
+func (qr *QR) qRowTo(i int, dst []float64) {
+	c := blas64.General{
+		Rows:   1,
+		Cols:   len(dst),
+		Stride: len(dst),
+		Data:   dst,
+	}
+	c.Data[i] = 1 // C is the i-th unit vector
+
+	// Construct Qi from the elementary reflectors: Qi = C * (H(1) H(2) ... H(nTau))
+	work := []float64{0}
+	lapack64.Ormqr(blas.Right, blas.NoTrans, qr.qr.mat, qr.tau, c, work, -1)
+	work = getFloat64s(int(work[0]), false)
+	lapack64.Ormqr(blas.Right, blas.NoTrans, qr.qr.mat, qr.tau, c, work, len(work))
+	putFloat64s(work)
+}
+
 // T performs an implicit transpose by returning the receiver inside a
 // Transpose.
 func (qr *QR) T() Matrix {
@@ -98,7 +136,9 @@ func (qr *QR) factorize(a Matrix, norm lapack.MatrixNorm) {
 	lapack64.Geqrf(qr.qr.mat, qr.tau, work, len(work))
 	putFloat64s(work)
 	qr.updateCond(norm)
-	qr.updateQ()
+	if qr.q != nil {
+		qr.q.Reset()
+	}
 }
 
 func (qr *QR) updateQ() {
@@ -149,7 +189,7 @@ func (qr *QR) RTo(dst *Dense) {
 		dst.ReuseAs(r, c)
 	} else {
 		r2, c2 := dst.Dims()
-		if c != r2 || c != c2 {
+		if r != r2 || c != c2 {
 			panic(ErrShape)
 		}
 	}
@@ -192,6 +232,10 @@ func (qr *QR) QTo(dst *Dense) {
 			panic(ErrShape)
 		}
 	}
+
+	if qr.q == nil || qr.q.IsEmpty() {
+		qr.updateQ()
+	}
 	dst.Copy(qr.q)
 }
 


Test patch:
diff --git a/mat/qr_test.go b/mat/qr_test.go
index b71bee56b..e9ec26aa2 100644
--- a/mat/qr_test.go
+++ b/mat/qr_test.go
@@ -18,9 +18,11 @@ func TestQR(t *testing.T) {
 	rnd := rand.New(rand.NewSource(1))
 	for _, test := range []struct {
 		m, n int
+		big  bool
 	}{
-		{5, 5},
-		{10, 5},
+		{m: 5, n: 5},
+		{m: 10, n: 5},
+		{m: 1e5, n: 3, big: true}, // Test that very tall matrices do not OoM.
 	} {
 		m := test.m
 		n := test.n
@@ -35,6 +37,13 @@ func TestQR(t *testing.T) {
 
 		var qr QR
 		qr.Factorize(a)
+		if test.big {
+			_ = qr.At(0, 0)     // should not panic, even for big matrices
+			_ = qr.At(m-1, n-1) // should not panic, even for big matrices
+			// We cannot proceed past here for big matrices.
+			continue
+		}
+
 		var q, r Dense
 		qr.QTo(&q)
 
@@ -56,6 +65,23 @@ func TestQR(t *testing.T) {
 		if !EqualApprox(&got, &want, 1e-12) {
 			t.Errorf("QR does not equal original matrix. \nWant: %v\nGot: %v", want, got)
 		}
+
+		// Verify indirect QR.At()
+		got.Reset()
+		got.ReuseAs(m, n)
+		qr.q.Reset() // reset q matrix to force lazy computation
+		for i := 0; i < m; i++ {
+			for j := 0; j < n; j++ {
+				got.set(i, j, qr.At(i, j))
+			}
+		}
+
+		if !EqualApprox(a, &got, 1e-14) {
+			t.Errorf("m=%d,n=%d: A and QR (computed with QR.At()) are not equal", m, n)
+		}
+		if !EqualApprox(a.T(), got.T(), 1e-14) {
+			t.Errorf("m=%d,n=%d: Aᵀ and (QR)ᵀ (computed with QR.At()) are not equal", m, n)
+		}
 	}
 }