LapackException(1) while svd but not svdvals

General Usage
1

I have a matrix such that

julia> svdvals(AA)
1000-element Array{Float64,1}:
 31.962791231191257  
 31.607090491558335  
 31.255070632248668  
 22.374974612534174  
  ⋮                  
 22.342788128226456  
 22.328813384016993  
  0.01373845014731832

but svd(AA) fails with the exception

julia> svd(AA)
ERROR: LAPACKException(1)
Stacktrace:
 [1] chklapackerror at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.2/LinearAlgebra/src/lapack.jl:38 [inlined]
 [2] gesdd!(::Char, ::Array{Float64,2}) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.2/LinearAlgebra/src/lapack.jl:1602
 [3] #svd!#73(::Bool, ::typeof(svd!), ::Array{Float64,2}) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.2/LinearAlgebra/src/svd.jl:63
 [4] #svd#74 at ./none:0 [inlined]
 [5] svd(::Array{Float64,2}) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.2/LinearAlgebra/src/svd.jl:106
 [6] top-level scope at none:0

Is it possible to tweak svd to use a more robust algorithm for this seemingly difficult problem?

Background: AA is a matrix of mostly Fourier basis functions that are orthogonal to each other, with the exception of one column which is a signal to be decomposed by the Fourier basis in a total least-squares fashion.

Code to generate AA is below

T = 1000
t = 0:T-1
f = LinRange(0,0.5-1/length(t)/2,length(t)÷2)
y = sin.(t)
function check_freq(f)
    zerofreq = findfirst(iszero, f)
    zerofreq !== nothing && zerofreq != 1 && throw(ArgumentError("If zero frequency is included it must be the first frequency"))
    zerofreq
end
function get_fourier_regressor(t,f)
    zerofreq = check_freq(f)
    N  = length(t)
    Nf = length(f)
    Nreg = zerofreq === nothing ? 2Nf : 2Nf-1
    N >= Nreg || throw(ArgumentError("Too many frequency components $Nreg > $N"))
    A  = zeros(N,Nreg)
    sinoffset = Nf
    for fn=1:Nf
        if fn == zerofreq
            sinoffset = Nf-1
        end
        for n = 1:N
            phi        = 2π*f[fn]*t[n]
            A[n,fn]    = cos(phi)
            if fn != zerofreq
                A[n,fn+sinoffset] = -sin(phi)
            end
        end
    end
    A, zerofreq
end
A,z = get_fourier_regressor(t,f)
AA  = [A y]
2
1 Like
3

Thanks for the pointer, I can manually call gesvd! for now and it works well :slight_smile:

1 Like
4

That looks like a very badly conditioned problem. Is it possible to reformulate it? I don’t know where the problem comes from but if I understand correctly, AA'*AA should be mostly identity (from Fourier orthogonality), except for the last column/line. This means that you can invert it semi-analytically using a Schur complement, and from there find the dominant eigenvalue of inv(AA'AA).

1 Like
5

You are right, a large part of AA'AA is almost orthogonal in this case since I chose f in a particular way. I do however not want to force this choice of f. The application is least-squares spectral analysis, which if f is chosen such that all basis functions are exactly orthogonal can be solved with the fourier transform. There are a number of reasons one may want to select the frequency vector differently though. Having said that, there might very well be a way of preconditioning the problem like you mention, I will think about it over night :slight_smile:
Thanks!

1 Like
6

well, having a non-orthogonal basis means that there’s always a price to pay down the line, so might as well orthogonalize it first thing :slight_smile:

2 Likes