DynareJulia / FastLapackInterface.jl

MIT License
32 stars 9 forks source link

Allow QR factorization with no columns. #39

Closed manuelbb-upb closed 8 months ago

manuelbb-upb commented 8 months ago

LinearAlgebra allows for QR decomposition of matrices with column number 0. In FastLapackInterface this is prohibited by an assertion in resize! Removing the assertion still errors, because LAPACK does not like nb==0. In the Standard Lib, there is an early return to don't even call geqrt, see this line We could easily do the same.

MichelJuillard commented 8 months ago

Yes, I will allow for 0-column matrix. Do you happen to know why Julia implementation of geqrt permit a 0-column matrix but not a 0-row matrix ? Shall we allow for both in FastLapackInterface ?

manuelbb-upb commented 8 months ago

Thanks! Regarding zero-row matrices in Julia's standard lib, I have no idea why they are not permitted. But I am not really familiar with the internals, to be honest, just traced down this particular issue. At first I thought that it could cause problems with QRCompactWY, but it seems to work:

import LinearAlgebra as LA
A = Matrix{Float64}(undef, 0, 10)
T = Matrix{Float64}(undef, 0, 0)
qr = LA.QRCompactWY(A, T)

and qr is as expected... I would be in favor of allowing it, but can also open an issue in the Julia repo to invite discussion by people with a bit more insight.

MichelJuillard commented 8 months ago

I will also allow for 0-row matrix and empty matrix. But, yes please, raise the issue in the Julia repo or on slack #linear-algebra

MichelJuillard commented 8 months ago

I have pushed branch https://gitlhub.com/DynareJulia/FastLapackInterface.jl/tree/zerocolumnrow, that allow for zero column and zero row matrix (or both). Could you please look at the changes and the additional tests and tell me if it is what you had in mind?

manuelbb-upb commented 8 months ago

Looks good to me! The issue came up, when I tried to optimize a routine, where QR factorization is used to determine a “sufficiently independent” set of interpolation/regression points for fully-linear surrogate models. If the initial interpolation set is empty, then there are no columns in the point matrix. Should work with those changes now.

MichelJuillard commented 8 months ago

Closed with commit https://github.com/DynareJulia/FastLapackInterface.jl/commit/7c05eb9e35d95820f64b76527cb3858409dc9ff0

MichelJuillard commented 7 months ago

The corresponding issue in Julia Linear algebra has been resolved https://github.com/JuliaLang/julia/issues/53451