1

通常の最小自乗回帰を解くための高速座標降下アルゴリズムを作成しようとしています。次ジュリア・コードは動作しますが、それはそんなにメモリに私は私がなぜジュリアはあまりにも多くのメモリを割り当てていますか?

BenchmarkTools.Trial: 
    memory estimate: 65.94 mb 
    allocs estimate: 151359 
    -------------- 
    minimum time:  19.316 ms (16.49% GC) 
    median time:  20.545 ms (16.60% GC) 
    mean time:  22.164 ms (16.24% GC) 
    maximum time:  42.114 ms (10.82% GC) 
    -------------- 
    samples:   226 
    evals/sample:  1 
    time tolerance: 5.00% 
    memory tolerance: 1.00% 

取得@benchmark OLS_cd(X, y)を使用して

n = 100 
p = 75 
σ = 0.1 
β_nz = float([i*(-1)^i for i in 1:10]) 

β = append!(β_nz,zeros(p-length(β_nz))) 
X = randn(n,p); X .-= mean(X,1); X ./= sqrt(sum(abs2(X),1)) 
y = X*β + σ*randn(n); y .-= mean(y); 

を生成するデータに

function OLS_cd{T<:Float64}(A::Array{T,2}, b::Array{T,1}, tolerance::T=1e-12) 
    N,P = size(A) 
    x = zeros(P) 
    r = copy(b) 
    d = ones(P) 
    while sum(d.*d) > tolerance 
     @inbounds for j = 1:P 
      d[j] = sum(A[:,j].*r) 
      x[j] += d[j] 
      r -= d[j]*A[:,j] 
     end 
    end 
    return(x) 
end 

を割り当てている理由私は理解していませんpが大きくなるにつれ、OLSの問題はより難しくなります。私はpを大きくして実行する必要があることに気付きました。

なぜ、それぞれwhileループを通過すると、より多くのメモリが割り当てられますか?私の目には、私の操作のすべてが適切であり、タイプが明確に指定されているようです。

プロファイリング中に何も飛び出さなかったが、有用であればその出力も投稿できる。

更新: としてベクトル化操作を使用することによって引き起こされる一時的なアレイが原因であった、以下に指摘。以下は、余分な割り当てを排除し、かなり迅速に実行します:

function OLS_cd_unrolled{T<:Float64}(A::Array{T,2}, b::Array{T,1}, tolerance::T=1e-12) 
    N,P = size(A) 
    x = zeros(P) 
    r = copy(b) 
    d = ones(P)  
    while norm(d,Inf) > tolerance 
     @inbounds for j = 1:P 
      d[j] = 0.0; @inbounds for i = 1:N d[j] += A[i,j]*r[i] end 
      @inbounds for i = 1:N r[i] -= d[j]*A[i,j] end 
      x[j] += d[j] 
     end 
    end 
    return(x) 
end 
+0

マイナーポイントは、Float64でなければならないと言ったときになぜパラメータ化を使用するのですか?ただそれをT <:Realにするのはもっと強力ではないでしょうか?また、ゼロ(P)と1(P)は正しいタイプを取りません、ゼロ(T、P)と1(T、P)でなければならず、内側のループでd [j] =ゼロ(T)。 –

+0

それは良い点です。私はまだJuliaの型指定のものでかなり荒いです。 – Rory

+0

また、配列からd [j]にアクセスする理由と、なぜ2つのループがあるのですか?どうしてですか? 'dj =ゼロ(T); @inbounds for i = 1:N; aij = A [i、j]; dj + = aij * r [i]である。 r [i] - = dj * aij;終わり ; x [j] + = dj;私は数学者ではないので、ここではすべて濡れているかもしれません。 –

答えて

4

A[:,j]は、コピーではなく、ビューを作成します。 @view A[:,j]またはview(A,:,j)を使用します。

r -= d[j]*A[:,j]r .= -.(r,d[j]*A[:.j])に分解すると、さらに一時的なものを取り除くことができます。 @ LutfullahTomakは、sum(A[:,j].*r)dot(view(A,:,j),r)としてデベクタケートしてそこにすべての一時的なものを取り除くべきだと述べた。中置演算子を使用するには、view(A,:,j)⋅rのように\cdotを使用できます。

コピーとビューについて、ベクトル化がテンポラリアレイをどのように引き起こすかについてお読みください。それは、ベクトル化された操作が発生すると、出力として新しいベクトルを作成する必要があるということです。代わりに、既存のベクトルに書きたいと思っています。 r = ...は配列を変更するので、配列を作成する式にはr = exが追加され、新しい配列が作成され、その配列にはrがポイントされます。 r .= exは、配列rの値を式の値で置き換えます。前者は一時的なものを割り当て、後者は一時的なものを割り当てます。このアイデアの繰り返し適用は、すべての一時的なものがどこから来るかです。

+0

私は '。='が0.6のためだと思っていましたが、それは0.5で私の驚きです。私は 'r =( - )。(r、(*)。(d [j]、view(A、:、5)))'を実行しました。私はそれが '( - )。(...)'であるべきだと思う。 –

+0

私はかっこが必要だとは思わない。しかし、スレッドのFYIとして、 '。+'の '.'が融合し、そのすべてが月曜日にマージされています。 [これは問題です](https://github.com/JuliaLang/julia/pull/17623#issuecomment-267784596)、数日間で夜行列車に乗ることになります。 Julia、私はそれを使用することをアドバイスしません) –

+0

私は仕事にビューアプローチを得ることができませんでしたが、forループにベクトル化された割り当てを展開するトリックをやったようだ。 – Rory

2

実際には、sum(d.*d)sum(A[:,j].*r)などは、インプレースではなく、一時的な配列を作る。..まず、sum(d.*d) == dot(d,d)と思うとsum(A[:,j].*r)は2つのテンポラリ配列になります。私は後者のためにdot(view(A,:,j),r)をやります。現在安定しているjulia(0.5)のバージョンはr -= d[j]*A[:,j]の短いバージョンではありませんので、ループを解く必要があります。

関連する問題