Well, if it’s correct to use a Cholesky factorization here it’s by far the most performant option given in this thread yet:

```
using Random
using RandomMatrices
function f1(A, Y, R)
(op1=A\Y,
op2=real.(sqrt(inv(A)))*R)
end
function f2(A::Hermitian, Y, R)
(op1 = A\Y,
op2 = sqrt(A)\R)
end
function f3(A::Hermitian, Y, R)
λ, ϕ = eigen(A)
(op1 = ϕ * Diagonal(1 ./ λ) * ϕ' * Y,
op2 = ϕ * Diagonal( 1 ./ sqrt.(λ)) * ϕ' * R)
end
function f4(A::Hermitian, Y, R)
C = cholesky(A)
(op1 = C \ Y,
op2 = C.L \ R)
end
```

```
let N = 100, Nx = 50
U = randn(Nx, N) #switched Nx and N so the code runs
A = I + U'*U
Y = randn(N, Nx) #switched Nx and N so the code runs
# Generate Random Matrices
d = Haar(1)
R = rand(d, N)
r1 = @btime f1($A, $Y, $R)
r2 = @btime f2($(Hermitian(A)), $Y, $R)
r3 = @btime f3($(Hermitian(A)), $Y, $R)
r4 = @btime f4($(Hermitian(A)), $Y, $R)
(same_results = Tuple(r1) .≈ Tuple(r4),
same_eigvalues = eigvals(r4.op2) ≈ eigvals(r1.op2),
same_singvalues = svd(r4.op2).S ≈ svd(r1.op2).S)
end
:RESULTS:
: 9.786 ms (31 allocations: 1.38 MiB)
: 1.780 ms (33 allocations: 705.03 KiB)
: 1.173 ms (31 allocations: 704.16 KiB)
: 117.878 μs (10 allocations: 273.80 KiB)
: (same_results = (true, false), same_eigvalues = false, same_singvalues = true)
:END:
```

As you can see, you don’t get the same answer for the second operation when you use the Cholesky L factor, but that doesn’t mean it isn’t equivalent. The resulting matrix appears to have different eigenvalues, but the same singular values as the other methods.