I’m implementing the Square-Root Kalman-Filter similar to this python implementation:
Python - Square-Root Kalman-Filter
The implementation calculates the QR-decomposition, but it only needs the upper triangular matrix R (see here). The Q is not needed. I searched the internet to see, if there was an efficient implementation to just calculate the matrix R and not Q and found this matlab page: Matlab qr
It states: To calculate the upper-triangular factor R
and permutation matrix P
, but avoid computing the orthogonal matrix Q
(which is often the most computationally expensive part of a call to qr
), you can specify B
as an empty matrix:
emptyB = zeros(size(S,1),0);
[~,R,P] = qr(S,emptyB);
I got curious and tried to find something similar in Julia. The closest I could get was the following:
A = randn(10,5)
B = zeros(10,1)
LAPACK.gels!('N', copy(A), B)
This results in the same output for R
as qr(A)
And indeed it is faster:
@btime LAPACK.gels!($'N', W, $B) setup=(W=copy(A)) evals=1 # 2.238 μs (7 allocations: 2.23 KiB)
@btime qr!(W) setup=(W=copy(A)) evals=1 # 3.772 μs (3 allocations: 608 bytes)
I have two questions:
- Why are there so many allocations even though
LAPACK.gels!
is an inplace function? - Is there a better way to calculate the matrix
R
? I tried to use a zero sizedB
as suggested by the matlab page:LAPACK.gels!('N', copy(A), zeros(10,0))
, but this does not provide the correct result.
I’d like to make the package publicly available as soon as it is finished.