Non reproducible result - Linear algebra

Specific Domains
,
1

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:

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:

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 is perfect.

Thank you very much

2

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

4 Likes
3

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 !

4

One other minor thing is that rather than doing

    ω = 2π*(freq .- freq[1]) # Reduced angular frequency
    X0 = exp.(-1im*ω*modelOrder'*Δt)

you can do

    ω = 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)

2 Likes
5

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.

6

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.

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, 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 , 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.

4 Likes
7

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.

8

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.

1 Like
9

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.