Support for non-square R matrices?

[njwfish]

I am wondering if it would be possible to add support for non-square R matrices for low rank A matrices.

As is the code of course allows this but it leads to huge numerical instabilities. It also doesn't "fail loudly" in the case where a singular column is being added. I only figured out this problem when I investigated some odd results I was getting (due to the numerical issues).

I am not sure it is possible to do Q-less updates in the low-rank case but if it is possible it'd be a really nice addition.

[njwfish]

I'll illustrate this and why it would be nice for rows then columns.

Columns

We can see that iteratively adding columns we get the same behaviour as qr:

A = rand(n, p)
init = 1
idx = 1:init
_, r = qr(A[:, idx])
for i in (init + 1):10
    idx = 1:i
    r = qraddcol(A[:, 1:(i - 1)], r, A[:, i])
    _, R = qr(A[:, idx])
    println(i, " ", norm(R' * R - r' * r))
end

> 2 6.153480596427405e-14
> 3 8.774555507236292e-14
> 4 1.2779896784284225e-13
> 5 1.5648080552598201e-13
> 6 1.8586499606283736e-13
> 7 2.1122966836905375e-13
> 8 2.5450922704198304e-13
> 9 2.756044971366047e-13
> 10 3.1800260559009874e-13

Adding singularity we get unpredictable behavior:

A = rand(n, p)
A[:, 3] = A[:, 4]
init = 1
idx = 1:init
_, r = qr(A[:, idx])
for i in (init + 1):10
    idx = 1:i
    r = qraddcol(A[:, 1:(i - 1)], r, A[:, i])
    _, R = qr(A[:, idx])
    println(i, " ", norm(R' * R - r' * r))
end

Rerunning this snippet sometimes we maintain precision, but other times we get LAPACKException(4) or, worst of all, we silently fail inducing large error in r' * r.

Rows

When n > p everything is fine with adding rows:

A = rand(100, 5)
init = 10
idx = 1:init
_, r = qr(A[idx, :])
for i in (init + 1):(init + 10)
    idx = 1:i
    r = qraddrow(r, A[i, :])
    _, R = qr(A[idx, :])
    println(i, " ", norm(R' * R - r' * r))
end
> 11 4.3737714791252564e-15
> 12 3.688882027804963e-15
> 13 5.9910909847502765e-15
> 14 3.9471512357794175e-15
> 15 5.4025784115714076e-15
> 16 4.905125646631375e-15
> 17 6.8510746500567485e-15
> 18 5.7902148899581444e-15
> 19 7.431033801933e-15
> 20 9.890336306082086e-15

But when p > n we have problems:

A = rand(100, 5)
init = 3
idx = 1:init
_, r = qr(A[idx, :])
for i in (init + 1):(init + 10)
    idx = 1:i
    r = qraddrow(r, A[i, :])
    _, R = qr(A[idx, :])
    println(i, " ", norm(R' * R - r' * r))
end
> 4 1.7222092200769785
> 5 2.7506485754003505
> 6 5.113888247250268
> 7 5.0845584882875015
> 8 7.376575144064207
> 9 8.304341679454847
> 10 9.177845742475348
> 11 10.290840634662828
> 12 10.533368921853537
> 13 11.04063738007243

Note that this never raises an error, it just silently fails while looking reasonable.

Conclusion

I just wanted to sketch out the fairly common patterns where the current code silently fails. I spent several hours over the last week debugging issues caused by these problems respectively. I think it would be great to either 1) somehow fail loudly or 2) update the code to run correctly for any size R matrix, rather than just square matrices.

I think this document lays out algorithms which could extend the methods in the current codebase. If thats correct I think it would be great to implement those algorithms here.