mat: calculate Q lazily when calling QR.ToQ

When a matrix is very tall, calculating Q will currently allocate a large Q at
the end of the factorisation, even if it is not going to be used. The eager
calculation was intended to prevent repeated re-calculation of Q when it is
used. So move the Q calculation to ToQ, but make it conditional on the stored
Q value being empty, and empty Q at the end of the factorisation.

mat/qr.go

@@ -98,23 +98,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()
-}
-
-func (qr *QR) updateQ() {
-	m, _ := qr.Dims()
-	if qr.q == nil {
-		qr.q = NewDense(m, m, nil)
-	} else {
-		qr.q.reuseAsNonZeroed(m, m)
+	if qr.q != nil {
+		qr.q.Reset()
 	}
-	// Construct Q from the elementary reflectors.
-	qr.q.Copy(qr.qr)
-	work := []float64{0}
-	lapack64.Orgqr(qr.q.mat, qr.tau, work, -1)
-	work = getFloat64s(int(work[0]), false)
-	lapack64.Orgqr(qr.q.mat, qr.tau, work, len(work))
-	putFloat64s(work)
 }
 
 // isValid returns whether the receiver contains a factorization.
@@ -192,9 +178,28 @@ func (qr *QR) QTo(dst *Dense) {
 			panic(ErrShape)
 		}
 	}
+	if qr.q == nil || qr.q.IsEmpty() {
+		qr.updateQ()
+	}
 	dst.Copy(qr.q)
 }
 
+func (qr *QR) updateQ() {
+	m, _ := qr.Dims()
+	if qr.q == nil {
+		qr.q = NewDense(m, m, nil)
+	} else {
+		qr.q.reuseAsNonZeroed(m, m)
+	}
+	// Construct Q from the elementary reflectors.
+	qr.q.Copy(qr.qr)
+	work := []float64{0}
+	lapack64.Orgqr(qr.q.mat, qr.tau, work, -1)
+	work = getFloat64s(int(work[0]), false)
+	lapack64.Orgqr(qr.q.mat, qr.tau, work, len(work))
+	putFloat64s(work)
+}
+
 // SolveTo finds a minimum-norm solution to a system of linear equations defined
 // by the matrices A and b, where A is an m×n matrix represented in its QR factorized
 // form. If A is singular or near-singular a Condition error is returned.

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,7 +37,15 @@ func TestQR(t *testing.T) {
 
 		var qr QR
 		qr.Factorize(a)
-		var q, r Dense
+
+		var r Dense
+		qr.RTo(&r)
+		if test.big {
+			// We cannot proceed past here for big matrices.
+			continue
+		}
+
+		var q Dense
 		qr.QTo(&q)
 
 		if !isOrthonormal(&q, 1e-10) {
@@ -49,8 +59,6 @@ func TestQR(t *testing.T) {
 			t.Errorf("m=%d,n=%d: Aᵀ and (QR)ᵀ are not equal", m, n)
 		}
 
-		qr.RTo(&r)
-
 		var got Dense
 		got.Mul(&q, &r)
 		if !EqualApprox(&got, &want, 1e-12) {