# Square-Root Kalman-Filter: Calculate only upper triangular of QR-decomposition
**URL:** https://discourse.julialang.org/t/square-root-kalman-filter-calculate-only-upper-triangular-of-qr-decomposition/55621
**Category:** Performance
**Tags:** kalman, qr
**Created:** [February 19, 2021, 2:59pm UTC](https://discourse.julialang.org/t/square-root-kalman-filter-calculate-only-upper-triangular-of-qr-decomposition/55621 "2021-02-19T14:59:45Z")
**Posts on this page:** 10
**Page:** 1
### Author:  [@zsoerenm](https://discourse.julialang.org/u/zsoerenm)
#### Post date: [February 19, 2021, 2:59pm UTC](https://discourse.julialang.org/t/square-root-kalman-filter-calculate-only-upper-triangular-of-qr-decomposition/55621/1 "2021-02-19T14:59:45Z")
I’m implementing the Square-Root Kalman-Filter similar to this python implementation:
[Python - Square-Root Kalman-Filter](https://github.com/rlabbe/filterpy/blob/master/filterpy/kalman/square_root.py)
The implementation calculates the QR-decomposition, but it only needs the upper triangular matrix R ([see here](https://github.com/rlabbe/filterpy/blob/master/filterpy/kalman/square_root.py#L210)). 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](https://de.mathworks.com/help/matlab/ref/qr.html)
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:
```julia
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:
```julia
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:
```julia
@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:
1. Why are there so many allocations even though `LAPACK.gels!` is an inplace function?
2. Is there a better way to calculate the matrix `R`? I tried to use a zero sized `B` 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.
---
### Author:  [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)
#### Post date: [February 19, 2021, 8:44pm UTC](https://discourse.julialang.org/t/square-root-kalman-filter-calculate-only-upper-triangular-of-qr-decomposition/55621/2 "2021-02-19T20:44:55Z")
Are you sure you actually need the square root form? I find that by simply symmetrizing the covariance matrix in the right places, you get almost the same numerical performance with a textbook implementation.
---
### Author:  [@BLI](https://discourse.julialang.org/u/BLI)
#### Post date: [February 19, 2021, 9:48pm UTC](https://discourse.julialang.org/t/square-root-kalman-filter-calculate-only-upper-triangular-of-qr-decomposition/55621/3 "2021-02-19T21:48:14Z")
I discussed the need of the square root KF with some industry people some 10 years ago (defence industry). Their point of view was that the square root KF is useful if one does calculations in Float32, but is not needed with Float64.
---
### Author:  [@zsoerenm](https://discourse.julialang.org/u/zsoerenm)
#### Post date: [February 20, 2021, 8:11am UTC](https://discourse.julialang.org/t/square-root-kalman-filter-calculate-only-upper-triangular-of-qr-decomposition/55621/4 "2021-02-20T08:11:22Z")
According to this paper [The square-root unscented Kalman filter for state and parameter-estimation](https://ieeexplore.ieee.org/document/940586) the square root form for the Unscented-Kalman-Filter has the following advantages:
- better numerical properties (compared to the UKF)
- guarantees positive semi-definiteness of the underlying state covariance
In my post above I asked for the standard Kalman-Filter, because I’d like to be able to switch from standard Kalman-Filter to the UKF or vice-versa depending on whether the process model or the measurement model is linear or non-linear.
Nonetheless, the SR-UKF also needs the matrix `R` from the `QR`-decomposition.
---
### Author:  [@BLI](https://discourse.julialang.org/u/BLI)
#### Post date: [February 20, 2021, 8:40am UTC](https://discourse.julialang.org/t/square-root-kalman-filter-calculate-only-upper-triangular-of-qr-decomposition/55621/5 "2021-02-20T08:40:47Z")
A classic book on factorization methods and the Kalman filter is Gerald J. Bierman’s Academic Press book from 1977; [Factorization Methods for Discrete Sequential Estimation | Gerald J. Bierman (Eds.) | download (1lib.sk)](https://1lib.sk/book/827556/53ee74) . This was, of course, before the UKF, so there are probably newer books, but perhaps Bierman’s book contains some motivation. He was at the Jet Propulsion Lab in Pasadena at the time.
---
### Author:  [@KalmanGuy](https://discourse.julialang.org/u/KalmanGuy)
#### Post date: [August 3, 2024, 3:17pm UTC](https://discourse.julialang.org/t/square-root-kalman-filter-calculate-only-upper-triangular-of-qr-decomposition/55621/6 "2024-08-03T15:17:24Z")
The UD algorithms from Bierman are the most robust, however, a Cholesky algorithm of the unscented filter is easiest to implement. Has anyone noticed the bug in the unscented filter equations? As documented, it will estimate observable parameters.
---
### Author:  [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)
#### Post date: [August 4, 2024, 5:24am UTC](https://discourse.julialang.org/t/square-root-kalman-filter-calculate-only-upper-triangular-of-qr-decomposition/55621/7 "2024-08-04T05:24:49Z")
> [@zsoerenm](#):
>
> I have two questions:
>
> 1. Why are there so many allocations even though `LAPACK.gels!` is an inplace function?
> 2. Is there a better way to calculate the matrix `R`? I tried to use a zero sized `B` as suggested by the matlab page: `LAPACK.gels!('N', copy(A), zeros(10,0))`, but this does not provide the correct result.
Have you found answers to your questions?
---
### Author:  [@zsoerenm](https://discourse.julialang.org/u/zsoerenm)
#### Post date: [August 4, 2024, 7:22am UTC](https://discourse.julialang.org/t/square-root-kalman-filter-calculate-only-upper-triangular-of-qr-decomposition/55621/8 "2024-08-04T07:22:23Z")
Yes, I implemented a non allocating `gels!` function here [KalmanFilters.jl/src/gels.jl at 1353334a50d9efc6b67c060557a0c8742a2da5b0 · JuliaGNSS/KalmanFilters.jl · GitHub](https://github.com/JuliaGNSS/KalmanFilters.jl/blob/1353334a50d9efc6b67c060557a0c8742a2da5b0/src/gels.jl#L11)
However, I haven’t found a way to use a zero sized `B`, but maybe that’s fine and wouldn’t have made a difference anyway.
---
### Author:  [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)
#### Post date: [August 7, 2024, 6:02pm UTC](https://discourse.julialang.org/t/square-root-kalman-filter-calculate-only-upper-triangular-of-qr-decomposition/55621/9 "2024-08-07T18:02:13Z")
Maybe it should be part of [FastLapackInterface.jl](https://github.com/DynareJulia/FastLapackInterface.jl).
@MichelJuillard, What do you think?
---
### Author:  [@MichelJuillard](https://discourse.julialang.org/u/MichelJuillard)
#### Post date: [August 8, 2024, 1:23pm UTC](https://discourse.julialang.org/t/square-root-kalman-filter-calculate-only-upper-triangular-of-qr-decomposition/55621/10 "2024-08-08T13:23:54Z")
Absolutely. Thanks for suggesting it. See [Add xGELS function · Issue #42 · DynareJulia/FastLapackInterface.jl · GitHub](https://github.com/DynareJulia/FastLapackInterface.jl/issues/42)
Thanks for suggesting it