# Non reproducible result - Linear algebra
**URL:** https://discourse.julialang.org/t/non-reproducible-result-linear-algebra/132986
**Category:** Specific Domains
**Tags:** question, linearalgebra
**Created:** [October 8, 2025, 4:28pm UTC](https://discourse.julialang.org/t/non-reproducible-result-linear-algebra/132986 "2025-10-08T16:28:24Z")
**Posts on this page:** 9
**Page:** 1
### Author:  [@Maucejo](https://discourse.julialang.org/u/Maucejo)
#### Post date: [October 8, 2025, 4:28pm UTC](https://discourse.julialang.org/t/non-reproducible-result-linear-algebra/132986/1 "2025-10-08T16:28:25Z")
Hi all,
I am implementing some algorithms for modal parameters extraction for the next release of `StructuralVibration.jl`. One of the method implemented is the Least-Squares Complex Frequency method (LSCF). The code often produce the right results, but for a reason I don’t understand it sometimes fails.
Here is a MWE to reproduce the result:
```julia
using StructuralVibration
# Structure parameters of the beam
L = 1. # Length
b = 0.03 # Width
h = 0.01 # Thickness
S = b*h # Cross-section area
Iz = b*h^3/12 # Moment of inertia
# Material parameters
E = 2.1e11 # Young's modulus
ρ = 7850. # Density
ξ = 0.01 # Damping ratio
# Mesh
xexc = 0:0.05:L
xm = xexc[2]
# Mode calculation - Simply supported boundary conditions
beam = Beam(L, S, Iz, E, ρ)
fmax = 500.
ωn, kn = modefreq(beam, 2fmax)
ϕexc = modeshape(beam, kn, xexc)
ϕm = modeshape(beam, kn, xm)
# FRF calculation
freq = 1.:0.1:fmax
prob = ModalFRFProblem(ωn, ξ, freq, ϕm, ϕexc)
H = solve(prob; ismat = true).u
p_lscf = lscf(H, freq, 20)
```
Here is the function `lscf`:
```julia
using LinearAlgebra, Polynomials
function lscf(frf, freq, order)
# FRF post-processing - Reshape FRF matrix
(nm, ne, nf) = size(frf)
FRF = reshape(frf, (nm*ne, nf))
# LSCF computation
ω = 2π*(freq .- freq[1]) # Reduced angular frequency
Δt = 1/(2(freq[end] - freq[1])) # Sampling interval
modelOrder = 0:order
nmodel = order + 1
# Reduced normal equation computation
M = similar(freq, nmodel, nmodel)
Rk = similar(M)
Sk = similar(M)
Tk = similar(M)
X0 = exp.(-1im*ω*modelOrder'*Δt)
Yk = similar(X0)
for Hk in eachrow(FRF)
# Build Yk
@. Yk = -Hk*X0
# Build Rk, Sk, Tk
Rk .= real(X0'X0)
Sk .= real(X0'Yk)
Tk .= real(Yk'Yk)
# Accumulate M
M .+= (Tk .- Sk'*(Rk\Sk))
end
# Solve the normal equation
A = -M[1:order, 1:order]
b = M[1:order, nmodel]
α = [A\b; 1.]
# Compute the poles
V = roots(Polynomial(α))
p = -log.(V)/Δt
# Poles that do not appear as pairs of complex conjugate numbers with a positive real part or that are purely real are suppressed. For complex conjugate poles, only the pole with a positive imaginary part is retained.
valid_poles = p[@. imag(p) < 0. && real(p) ≤ 0. && !isreal(p)]
poles = intersect(p, conj.(valid_poles))
# To avoid modifying the stability of the dynamic system in case of frequency offset, we set real(poles_new) = real(poles_old). This means that dn_new = dn_old * fn_old / fn_new.
if freq[1] != 0.
fn, ξn = poles2modal(poles)
poles .= modal2poles(fn .+ freq[1], @. ξn*fn/(fn + freq[1]))
end
sort!(poles, by = abs)
return poles
end
function modal2poles(fn, ξn)
ωn = 2π*fn
return @. -ξn*ωn + 1im*ωn*√(1 - ξn^2)
end
function poles2modal(poles)
ωn = abs.(poles)
fn = ωn/(2π)
ξn = @. -real(poles)/ωn
return fn, ξn
end
```
I suspect that my problem is related to the matrix inversions I have to do, but I am not sure why.
If you have any idea or suggestions to help me to make the code more robust numerically, I would really appreciate.
For checking the math behind, this [paper](https://www.researchgate.net/publication/265423092_A_poly-reference_implementation_of_the_least-squares_complex_frequency-domain_estimator) is perfect.
Thank you very much
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [October 8, 2025, 5:22pm UTC](https://discourse.julialang.org/t/non-reproducible-result-linear-algebra/132986/2 "2025-10-08T17:22:29Z")
You likely don’t want to be using the normal equations. They lose a lot of numeric stability compared to just solving directly (which will likely use a QR factorization).
Specifically, rather than using `Sk'*(Rk\Sk)` you can use `X0\Yk`
---
### Author:  [@Maucejo](https://discourse.julialang.org/u/Maucejo)
#### Post date: [October 8, 2025, 7:17pm UTC](https://discourse.julialang.org/t/non-reproducible-result-linear-algebra/132986/3 "2025-10-08T19:17:01Z")
Thank you @Oscar_Smith. I rewrite my function to solve the problem using either the SVD (and eliminating the smallest singular values below some threshold) or the QR factorization. In both cases, the numerical stability is far more greater !
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [October 8, 2025, 7:43pm UTC](https://discourse.julialang.org/t/non-reproducible-result-linear-algebra/132986/4 "2025-10-08T19:43:10Z")
One other minor thing is that rather than doing
```julia-auto
ω = 2π*(freq .- freq[1]) # Reduced angular frequency
X0 = exp.(-1im*ω*modelOrder'*Δt)
```
you can do
```julia-auto
ω = 2*(freq .- freq[1]) # Reduced angular frequency
X0 = cispi.(-ω*modelOrder'*Δt)
```
where `cispi(x)` computes `exp(im*pi*x)` (but does so more accurately and quicker)
---
### Author:  [@Maucejo](https://discourse.julialang.org/u/Maucejo)
#### Post date: [October 9, 2025, 5:19am UTC](https://discourse.julialang.org/t/non-reproducible-result-linear-algebra/132986/5 "2025-10-09T05:19:11Z")
I am confused because I am translating a Matlab code into Julia. The flow chart for the backslash operator is different in Matlab than in Julia. In my case, Matlab uses a QR factorization, while Julia uses an algorithm similar to `lsqminnorm`.
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [October 9, 2025, 5:54am UTC](https://discourse.julialang.org/t/non-reproducible-result-linear-algebra/132986/6 "2025-10-09T05:54:40Z")
> [@Maucejo](#):
>
> In my case, Matlab uses a QR factorization, while Julia uses an algorithm similar to `lsqminnorm`.
They both use a QR factorization. Julia’s `A \ b` also uses (pivoted) QR factorization of `A` for a “tall” (overdetermined) system A, and gives the least-square solution in this case, much like Matlab.
However, they differ for an underdetermined system (“wide” A), or an overdetermined A that is not full column rank. In both of these cases the least-square solution is not unique. In such cases, Julia returns the minimum-norm solution, whereas Matlab returns a “basic” solution IIRC. This boils down to using a QR factorization in different ways.
> [@Maucejo](#):
>
> I rewrite my function to solve the problem using either the SVD (and eliminating the smallest singular values below some threshold) or the QR factorization.
If you are going to explicitly regularize the problem by truncating the SVD, that is different from either of the above, unless your threshold is proportional to machine precision. If you just do `svd(A) \ b`, it will [already drop singular values near machine precision](https://github.com/JuliaLang/LinearAlgebra.jl/blob/db48f0dd6f6a33aee8068718e1338628a9d388d1/src/svd.jl#L267-L272), so you don’t need to manually truncate if that’s all you want. (This is not explicitly documented yet, but will be documented in the next release by the PR below.)
Note also that if you are going to use the SVD A = U \Sigma V^\*, it is important to apply the pseudo-inverse as A^+ b = V( \Sigma^+ (U^\* b)), with U^\* b computed first, rather than explicitly computing the `pinv` matrix. The latter is numerically unstable. See [stable pinv least-squares via rtol,atol in svd, ldiv! by stevengj · Pull Request #1387 · JuliaLang/LinearAlgebra.jl · GitHub](https://github.com/JuliaLang/LinearAlgebra.jl/pull/1387) , which will also give a convenient syntax (once it appears in a Julia release) for putting a threshold on the SVD: it introduces new syntaxes like `SVD(A; rtol, atol) \ b` for conveniently applying a truncated SVD.
Also, using a truncated SVD gives the minimum-norm solution in cases where the solution is not unique, just like `A \ b`.
---
### Author:  [@Maucejo](https://discourse.julialang.org/u/Maucejo)
#### Post date: [October 9, 2025, 10:11am UTC](https://discourse.julialang.org/t/non-reproducible-result-linear-algebra/132986/7 "2025-10-09T10:11:08Z")
I implemented the `lscf` function using the naive implementation presented in this thread in Python and MATLAB, and I didn’t observe the numerical instability issue. I’m curious what kind of “wizardry” is used in the other languages.
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [October 9, 2025, 10:42am UTC](https://discourse.julialang.org/t/non-reproducible-result-linear-algebra/132986/8 "2025-10-09T10:42:27Z")
Matlab will automatically (and silently) switch to a least-square solution if you do `S \ b` with a singular _square_ matrix, whereas Julia will give an error. Not sure what Python function you were using.
However, they are probably losing accuracy. If you solve a least-squares problem \min\_x \Vert Ax -b \Vert by the normal-equations approach (A^T A)^{-1} (A^T b), you are sunk (you’ve squared the condition number, i.e. doubled the number of digits lost to roundoff error) as soon as you compute A^T A. Nothing in these languages will save you from that.
---
### Author:  [@Maucejo](https://discourse.julialang.org/u/Maucejo)
#### Post date: [October 9, 2025, 10:52am UTC](https://discourse.julialang.org/t/non-reproducible-result-linear-algebra/132986/9 "2025-10-09T10:52:31Z")
Thank you for this explanation. In Python, I use the `solve` function from scipy.linalg.
I prefer the julia approach. I want to receive an error when needed to avoid hidden numerical problems.