# LapackException(1) while svd but not svdvals
**URL:** https://discourse.julialang.org/t/lapackexception-1-while-svd-but-not-svdvals/23787
**Category:** General Usage
**Tags:** linearalgebra
**Created:** [May 2, 2019, 5:57pm UTC](https://discourse.julialang.org/t/lapackexception-1-while-svd-but-not-svdvals/23787 "2019-05-02T17:57:47Z")
**Posts on this page:** 6
**Page:** 1
### Author:  [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)
#### Post date: [May 2, 2019, 5:57pm UTC](https://discourse.julialang.org/t/lapackexception-1-while-svd-but-not-svdvals/23787/1 "2019-05-02T17:57:47Z")
I have a matrix such that
```julia
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
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
```julia
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]
```
---
### Author:  [@carstenbauer](https://discourse.julialang.org/u/carstenbauer)
#### Post date: [May 2, 2019, 6:28pm UTC](https://discourse.julialang.org/t/lapackexception-1-while-svd-but-not-svdvals/23787/2 "2019-05-02T18:28:11Z")
> [@baggepinnen](#):
>
> Is it possible to tweak `svd` to use a more robust algorithm for this seemingly difficult problem?
>
>
> This PR adds an \`alg\` keyword to \`svd\` and \`svd!\` which allows the user to switc…h between a divide-and-conquer algorithm (default, \`LAPACK.gesdd!\`) and a simple algorithm (\`LAPACK.gesvd!\`).
>
> In the future, \[Jacobi methods\](https://discourse.julialang.org/t/calling-lapacks-jacobi-svd-dgesvj/20602) could also be added. (#31080)
>
> See #31023 and https://discourse.julialang.org/t/svd-better-default-to-gesvd-instead-of-gesdd/20603 for related discussions.
>
> Closes #31023.
---
### Author:  [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)
#### Post date: [May 2, 2019, 6:45pm UTC](https://discourse.julialang.org/t/lapackexception-1-while-svd-but-not-svdvals/23787/3 "2019-05-02T18:45:52Z")
Thanks for the pointer, I can manually call `gesvd!` for now and it works well 🙂
---
### Author:  [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)
#### Post date: [May 2, 2019, 6:49pm UTC](https://discourse.julialang.org/t/lapackexception-1-while-svd-but-not-svdvals/23787/4 "2019-05-02T18:49:38Z")
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)`.
---
### Author:  [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)
#### Post date: [May 2, 2019, 6:56pm UTC](https://discourse.julialang.org/t/lapackexception-1-while-svd-but-not-svdvals/23787/5 "2019-05-02T18:56:50Z")
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 🙂
Thanks!
---
### Author:  [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)
#### Post date: [May 2, 2019, 6:59pm UTC](https://discourse.julialang.org/t/lapackexception-1-while-svd-but-not-svdvals/23787/6 "2019-05-02T18:59:04Z")
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 🙂