Hi everyone, I am trying to improve my codeβs efficiency. I locate a line that might be too expensive, but I failed to figure out how to improve it.

In brief, I need to save the inverse of an N by N covariance matrix of N realizations of a Gaussian process with Matern covariogram. In the following code chunk, I define a function `Maternlu`

to compute the Matern covariance between two points separated by distance x. `coords`

records the coordinates of all realizations and matrix `M`

stores what I need.

```
using Distributions
using LinearAlgebra
using Distances
using SpecialFunctions
function Maternlu(x; Ξ½=1.0, Ο = 6.0, Ο2 = 1.0)
# Matern function
if x == 0.0
Ο2
else
a = Ο * x;
Ο2 * 2^(1 - Ξ½) / gamma(Ξ½) * (a)^Ξ½ * besselk(Ξ½, a)
end
end
N = 400;
coords = rand(2, N); #coordinates of observations
M = inv(cholesky(Symmetric(Maternlu.(
UpperTriangular(pairwise(Euclidean(), coords, dims = 2)),
Ξ½ = 1.0, Ο = 6.0))));
```

The last line requires around 6MB but the matrix `M`

itself only requires 1.2MB. Moreover, I only need the upper-triangular part of `M`

.

```
@benchmark M = inv(cholesky(Symmetric(Maternlu.(
UpperTriangular(pairwise(Euclidean(), coords, dims = 2)),
Ξ½ = 1.0, Ο = 6.0))))
BenchmarkTools.Trial: 63 samples with 1 evaluation.
Range (min β¦ max): 39.158 ms β¦ 83.481 ms β GC (min β¦ max): 0.00% β¦ 3.10%
Time (median): 80.759 ms β GC (median): 0.00%
Time (mean Β± Ο): 80.414 ms Β± 5.400 ms β GC (mean Β± Ο): 0.27% Β± 0.73%
βββ
βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ β
39.2 ms Histogram: frequency by time 83.3 ms <
Memory estimate: 6.10 MiB, allocs estimate: 79819.
```

I think the ideal code will be like

- preallocate an N by N upper-triangular matrix
- overwrite it with the upper-triangular part of the distance matrix among the N locations
- overwrite it with the upper-triangular part of the Matern covariance
- use function like
`cholesky!`

to overwrite this matrix with the Cholesky decomposition - overwrite it with the upper-triangular part of the inverse of the Matern covariance

But I donβt know how to preallocate an upper-triangular matrix and it seems to me that the current `cholesky!`

doesnβt support an upper-triangular matrix input. Any thoughts about what I can do to improve my code? Many thanks!!!