Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

mat: calculate Q lazily when calling QR.ToQ #1970

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

mat: calculate Q lazily when calling QR.ToQ #1970

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
37 changes: 21 additions & 16 deletions mat/qr.go
View file
Open in desktop
Expand Up @@ -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()
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm afraid this doesn't work because q is currently needed in the At method. The test passes just because QTo is called in the test before EqualApprox (which calls At).

We could reconstruct the necessary row of Q in each call to At which is terrible but hopefully nobody uses QR as Matrix except in the call to Solve?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The other alternative is to lazily calculate Q for At the same way it is for ToQ, with a warning that it may cause OoM. I think your approach is probably better. A warning in the docs that it will be extremely inefficient should be enough.

Copy link

@tvkn tvkn Apr 18, 2024 •
edited

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We could lazily compute q in the At() method if q is nil or empty, in the same way it is done in the Qto.

But I guess I'm also a bit confused as to why we need q at all in At() in the first place. When we factorize a matrix a, qr stores a copy of that matrix. Why not return directly qr.At()?

Isn't https://github.com/gonum/gonum/blame/2ad11cabb395b96efc5b67fa1b64480762d9e703/mat/qr.go#L46 also faulty?
I'm expected we should return Q*R at element (i,j), but instead it seems we're returning Q*A at element (i,j), with A = Q*R, the matrix we factorize.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The other alternative is to lazily calculate Q for At the same way it is for ToQ.

Agreed. If one element is needed, then most likely all are.

qr stores a copy of that matrix. Why not return directly qr.At()?

That copy is overwritten by lapack64.Geqrf which efficiently stores both Q and R in the same matrix.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Another alternative is to indeed store a in Factorize and compute Q in QTo without storing it.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think we would need to copy a in this situation.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In the end I think it's best to lazily compute Q in At. QR.At should never be called in practice.

Copy link
Member Author

@kortschak kortschak May 20, 2024 •
edited

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm thinking that @tvkn's suggestion of conditionally computing element of QR for At via the relevant row of Q if Q is not already calculated, otherwise falling through to the code that is there now.

diff --git a/mat/qr.go b/mat/qr.go
index d4787bc3..7caa2f20 100644
--- a/mat/qr.go
+++ b/mat/qr.go
@@ -41,6 +41,11 @@ func (qr *QR) At(i, j int) float64 {
 		panic(ErrColAccess)
 	}
 
+	if qr.q == nil || qr.q.IsEmpty() {
+		// Calculate Q's row i;
+		// Return Q[i,]·R[,j].
+	}
+
 	var val float64
 	for k := 0; k <= j; k++ {
 		val += qr.q.at(i, k) * qr.qr.at(k, j)

This is safer than just lazily calculating Q which has the risk of someone thinking the code is fine until a tall matrix is used down the track. At worst with this approach there will be a perf hit, but not a crash.

I don't have time to do this, so if someone wants to pick this up. I'm happy for that to happen.

}

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.
Expand Down Expand Up @@ -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.
Expand Down
18 changes: 13 additions & 5 deletions mat/qr_test.go
View file
Open in desktop
Expand Up @@ -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
Expand All @@ -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) {
Expand All @@ -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) {
Expand Down