Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

L\b with cholesky L returning stackoverflow error #16196

Closed
cdipaolo opened this issue May 4, 2016 · 7 comments
Closed

L\b with cholesky L returning stackoverflow error #16196

cdipaolo opened this issue May 4, 2016 · 7 comments

Comments

@cdipaolo
Copy link

cdipaolo commented May 4, 2016

Reference #3630.

I have a matrix K_ with dimensions 252x252. This matrix is positive semidefinite (a gram matrix from a mercer kernel). For numerical purposes I am trying not to use

> pinv(K_)*tsla

252-element DataArrays.DataArray{Float64,1}:
 0.289256
 0.283591
 0.265911
 0.284096
 0.281966
 0.286117
 0.295026
 0.290874
 0.28921 
 0.291895
 0.293412
 0.29149 
 0.286033
        
 0.296586
 0.298044
 0.28968 
 0.293612
 0.294002
 0.301217
 0.300779
 0.302519
 0.299725
 0.298498
 0.293113
 0.288902

which does work. Because K_ is positive semidefinite I am trying to compute the Cholesky factorization K_ = L*L' first and then solve with L'\(L\tsla). This, however, returns a stack overflow error:

> L = chol(K_)
> L'\(L\tsla)

LoadError: StackOverflowError:
while loading In[92], in expression starting on line 2

What is interesting is that this is not happening in general. If I do, with bigger dimensions,

> A = randn(500,500)
> A = A'*A
> L = chol(A)
> b = randn(500,1)
> L\(L'\b)  pinv(A)*b

true

we have the correct result.

I am attaching my specific arrays here, created with writecsv()
matrices.zip

Is there a reason why this should happen?

@yuyichao
Copy link
Contributor

yuyichao commented May 4, 2016

What's your versioninfo()? Maybe not the same but sounds similar to JuliaLang/LinearAlgebra.jl#283

@cdipaolo
Copy link
Author

cdipaolo commented May 4, 2016

> versioninfo()

Julia Version 0.4.5
Commit 2ac304d* (2016-03-18 00:58 UTC)
Platform Info:
  System: Darwin (x86_64-apple-darwin14.5.0)
  CPU: Intel(R) Core(TM) i5-4260U CPU @ 1.40GHz
  WORD_SIZE: 64
  BLAS: libopenblas (DYNAMIC_ARCH NO_AFFINITY Haswell)
  LAPACK: libopenblas
  LIBM: libopenlibm
  LLVM: libLLVM-3.3

@yuyichao
Copy link
Contributor

yuyichao commented May 4, 2016

Hmmm, I can't reproduce on 0.4 linux...

               _
   _       _ _(_)_     |  A fresh approach to technical computing
  (_)     | (_) (_)    |  Documentation: http://docs.julialang.org
   _ _   _| |_  __ _   |  Type "?help" for help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 0.4.6-pre+18 (2016-04-01 03:37 UTC)
 _/ |\__'_|_|_|\__'_|  |  Commit 198644a* (33 days old release-0.4)
|__/                   |  x86_64-unknown-linux-gnu

julia> K_ = readcsv("K_.csv");

julia> tsla = readcsv("tsla.csv");

julia> L = chol(K_);

julia> L'\(L\tsla);

julia> 

Not on master either....

@cdipaolo
Copy link
Author

cdipaolo commented May 4, 2016

Oh that's interesting. Yeah I just reloaded and it works now, but if I construct the matrices directly it doesn't work.

> K_ = readcsv("../data/processed/matrices/K_.csv")
> tsla = readcsv("../data/processed/matrices/tsla.csv")
> L = chol(K_)
> L\(L'\tsla)  pinv(K_)*tsla
> t    = tsla_raw[:Date]
> tsla = tsla_raw[:Open]
> t_   = collect(minimum(t):0.25:maximum(t))
> println("==> Run GP Regression over TSLA Open Prices")
> println("--  between May 5, 2015 and May 4, 2016")
> σ = 2.
> l = 10.
> se = SEKernel(l,σ)
> σ_y = 28.
> println("==> Construct Kernels")
> K_ = K(se, t, t)
> K_ = K_ + σ_y^2*eye(K_)
> k_ = K(se, t, t_)
> dk_ = d_K(se,t,t_)
> k__ = K(se, t_,t_)
> println("--  done.")

> L = chol(K_)
> println(size(L))
> α = L\(L'\tsla)

==> Run GP Regression over TSLA Open Prices
--  between May 5, 2015 and May 4, 2016
==> Construct Kernels
--  done.
(252,252)
LoadError: StackOverflowError:
while loading In[132], in expression starting on line 26

@cdipaolo
Copy link
Author

cdipaolo commented May 4, 2016

OH this is weird. If I construct the matrices themselves and just load the tsla array from the csv file it works, so it seems as if the problem is with the b matrix in Ax=b itself.

t    = tsla_raw[:Date]
#tsla = tsla_raw[:Open]
tsla = readcsv("../data/processed/matrices/tsla.csv")
t_   = collect(minimum(t):0.25:maximum(t))

println("==> Run GP Regression over TSLA Open Prices")
println("--  between May 5, 2015 and May 4, 2016")

σ = 2.
l = 10.
se = SEKernel(l,σ)
σ_y = 28.

println("==> Construct Kernels")
K_ = K(se, t, t)
K_ = K_ + σ_y^2*eye(K_)

k_ = K(se, t, t_)
dk_ = d_K(se,t,t_)
d2k_ = d2_K(se,t,t_)
k__ = K(se, t_,t_)

println("--  done.")

L = chol(K_)
println(size(L))
α = L\(L'\tsla)

-----------------------------------------------------------------

==> Run GP Regression over TSLA Open Prices
--  between May 5, 2015 and May 4, 2016
==> Construct Kernels
--  done.
(252,252)
Out[133]:
252x1 Array{Float64,2}:
 0.289256
 0.283591
 0.265911
 0.284096
 0.281966
 0.286117
 0.295026
 0.290874
 0.28921 
 0.291895
 0.293412
 0.29149 
 0.286033
        
 0.296586
 0.298044
 0.28968 
 0.293612
 0.294002
 0.301217
 0.300779
 0.302519
 0.299725
 0.298498
 0.293113
 0.288902

@cdipaolo
Copy link
Author

cdipaolo commented May 4, 2016

The dimensions when I load the array from csv is 2 now, instead of 1:

> tsla = readcsv("../data/processed/matrices/tsla.csv")
Out[134]:
252x1 Array{Float64,2}:
 237.76
 234.1 
 221.0 
 235.99
 236.29
 240.11
 247.61
 244.82
 243.93
 247.0 
 248.43
 247.13
 243.03
      
 252.23
 253.12
 246.26
 248.99
 248.89
 253.01
 252.05
 252.75
 249.85
 248.14
 241.5 
 237.36
In [135]:

> tsla = tsla_raw[:Open]
Out[135]:
252-element DataArrays.DataArray{Float64,1}:
 237.76
 234.1 
 221.0 
 235.99
 236.29
 240.11
 247.61
 244.82
 243.93
 247.0 
 248.43
 247.13
 243.03
      
 252.23
 253.12
 246.26
 248.99
 248.89
 253.01
 252.05
 252.75
 249.85
 248.14
 241.5 
 237.36

@andreasnoack
Copy link
Member

Minimal example:

julia> using DataArrays

julia> UpperTriangular(randn(3,3))\DataArray(randn(3))
ERROR: StackOverflowError:
 in \ at linalg/generic.jl:323
 in \ at linalg/generic.jl:326 (repeats 65491 times)

but also

julia> UpperTriangular(randn(3,3))\sub(randn(3), [1,2,3])
ERROR: StackOverflowError:
 in \ at linalg/generic.jl:326 (repeats 74846 times)

We might want to change the signature in

function ($f){TA,TB,S}(A::Union{UnitUpperTriangular{TA,S},UnitLowerTriangular{TA,S}}, B::StridedVecOrMat{TB})
to AbstractVecOrMat.

More shortterm you might want to make sure that you convert the DataVector to a Vector. It will also be much faster.

Finally, notice that L is actually upper triangular.

andreasnoack added a commit that referenced this issue May 5, 2016
andreasnoack added a commit that referenced this issue May 17, 2016
andreasnoack added a commit that referenced this issue May 18, 2016
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

3 participants