トップ 追記

日々の流転


2014-04-30 [長年日記]

λ. IOモナドでヒープ割り当てを強制される件

Haskellでプログラムを書いていて、数年以上ずっと悩んでいる点として、

sumI :: UArray Int Int -> Int -> Int -> Int -> Int
sumI a i end ret
  | i <= end  = sumI a (i+1) end (ret + (a ! i))
  | otherwise = ret

みたいな関数だと、UArray Int Int -> Int# -> Int# -> Int# -> Int# という型のワーカ関数が生成されるので、結果を即消費するような場合にはヒープ割り当てが発生しないのに対して、

sumM :: IOUArray Int Int -> Int -> Int -> Int -> IO Int
sumM a i end !ret
  | i <= end = do
      x <- readArray a i
      sumM a (i+1) end (ret+x)
  | otherwise =
      return ret

みたいな関数だと、IOUArray Int Int -> Int# -> Int# -> Int# -> State# RealWorld -> (# State# RealWorld, Int# #) という型のワーカ関数ではなく、IOUArray Int Int -> Int# -> Int# -> Int# -> State# RealWorld -> (# State# RealWorld, Int #) という型のワーカーしか生成されず、Intが必ずヒープ割り当てされてしまう。

通常は大して問題にならないのだけど、凄まじい回数呼ばれる関数だと、それが積もり積もってプログラム全体のメモリ割り当ての数割を占めてしまったりして、そういう場合にはGCプレッシャーも馬鹿にならないので、出来ることなら1バイトたりともヒープ割り当てしたくないのである。 望む worker-wrapper transformation 結果を手で書いてしまえば良いのだけれど、GHC拡張に依存したコードを直書きするのも気が引けることもあって……

Tags: haskell
本日のツッコミ(全2件) [ツッコミを入れる]

ψ mkotha [ご存知かもしれませんがJoachim Breitnerさんがこの問題に取り組んでますね https://ghc.ha..]

ψ さかい [おお。知らなかったです。ありがとうございます!]


2013-12-31 [長年日記]

λ. 2013年振り返り

2014年を迎えるにあたって、2013年を簡単に振り返ってみる。

月毎の主なイベント

1月
Coursera での学習を開始。昨年末に買ったNexus7で初めてのAndroid経験&初代iPadから乗り換え。
2月
スキーで転けて靭帯損傷。
3月
『型システム入門 プログラミング言語と型の理論』出版と出版記念トークイベントMax-SAT Evaluation 2013 に参加。
4月
特になし
5月
Haskellで計算機代数勉強会を開催。Google Code Jam 2013 は去年に引き続いて今年も体調的に微妙で微妙。
6月
某自社の株主総会見てきた。
7月
スカイツリー昇った。ICSE'13勉強会でT芝チームとして発表。 奥歯抜歯Coursera MeetUp Tokyo に参加。
8月
ICFP Contest 2013 に例によって Team Sampou で参加。引越しと一人暮らし開始。SPLC2013の併設ワークショップFMSPLE 2013で発表。
9月
ウサギと同居開始。
10月
結婚(といっても婚姻届を出しただけで、結婚式はまだ)
11月
電機連合第34回技術者フォーラムに参加。2013川崎国際多摩川マラソンに参加。
12月
川崎グランシティモール実証実験開始

Coursera

前から少し気になっていたCourseraを1月に試してみて、最終的に以下の9コースを履修。

  1. Computing for Data Analysis by Roger D. Peng @ Johns Hopkins University
  2. Data Analysis by Jeff Leek @ Johns Hopkins University
  3. Linear and Discrete Optimization by Friedrich Eisenbrand @ École Polytechnique Fédérale de Lausanne
  4. Machine Learning by Andrew Ng @ Stanford University
  5. Developing Innovative Ideas for New Companies: The First Step in Entrepreneurship by James V. Green @ University of Maryland, College Park
  6. Model Thinking by Scott E. Page @ University of Michigan
  7. Maps and the Geospatial Revolution by Anthony C. Robinson @ The Pennsylvania State University
  8. Introduction to Computational Finance and Financial Econometrics by Eric Zivot @ University of Washington
  9. Introduction to Recommender Systems by Joseph A Konstan and Michael D Ekstrand @ University of Minnesota

自分のスキルのポートフォリオに追加したいと考えたものを中心に、それなりに戦略的に学習してみたつもり。上手くいった部分もあればそうでない部分もあるけれど。自分はデータ分析や統計に関しては大学1年のときに、ほぼ必修の位置づけだった「データ分析」の授業を途中で飽きて切ってしまってから、ほとんど学習しないままになってしまっていたので、その辺りを補完できたのは良かった。

そして、7月にCoursera MeetUp Tokyoに参加したところ、その縁で、メディアで取り上げらることになったのは、なかなか貴重な体験だった。

新しいプログラミング言語、環境

「毎年少なくとも一つの言語を学習する」という話があるけれど、今年は Coursera のコースで割と色々な言語を学べた。 Computing for Data AnalysisData AnalysisIntroduction to Recommender Systems では R を、Linear and Discrete Optimization ではPython (とSymPy)を、Machine Learningでは MATLAB/Octave に触れることができた。その後、そのうちのRとPythonは仕事でもそれなりに使うようになったので、結構役立っている。

それから、言語というほどではないけれど、Model Thinking では NetLogo に、Maps and the Geospatial Revolution では ArcGIS Online に触った。 あと、Introduction to Recommender Systems で使った LensKit は、アノテーションなどを活用したJavaの今風なDIを使ってて、ちょっとおぉと思ったりなんかもした。

『型システム入門 プログラミング言語と型の理論』出版

Benjamin C. Pierce の Types and Programming Languages (TAPL) の訳書『型システム入門 プログラミング言語と型の理論』を出版。前回のAlloy本こと『抽象によるソフトウェア設計−Alloyではじめる形式手法』が終わり、2011年9月からTAPL翻訳が本格始動したのだけれど、それから1年半ほどかかって、ついに出版にこぎ着いたのがこの3月で、非常に感慨深かった。また、共訳者、監訳者、編集者、レビュワー、皆素晴らしい人達ばかりで一緒に仕事をできて光栄だった。出版記念トークイベントという珍しい体験もしてしまったし。

靭帯損傷

2月にスキーで転けて靭帯損傷。内側側副靭帯損傷の2度で、MRIを初体験したり、しばらくギプス(正確にはギプスシーネ)と松葉杖で不便な生活をしたり、さらに結構な間リハビリをしたりとか。 ギブスしてると、立ってても座っててもバランスが悪くて、体中凝って仕方なかったのが辛かった。 後、しばらくギプスとサポーターをしてるだけで、関節って本当に曲がらなくなるんだなぁ。リハビリしてだんだん動くようになって、人間の関節や筋肉って本当によくできてるなぁと思った。 大変な体験だったけれど、周囲からは非常によくしてくれたし、結構学ぶことのあった体験だったと思う。

ランニング

ランニングは500kmの目標を立てていたけど、途中怪我で数ヶ月走れなかったこともあって、320kmほどに留まった。でも、我ながらよく続いたなぁと思うし、自分を誉めてあげたい。 また、大会に参加するという目標は2013川崎国際多摩川マラソンに参加するという形で達成。 目標がハーフマラソンだったのに対して8kmになってしまったので、ハーフマラソンを走るのは2014年の目標に持ち越し。

あと、昨年末にちゃんとしたランニングシューズを買ったけれど、今年は心拍計 Wahoo Blue HR とか、活動量計 Misfit Shine とか、ランニング用のスポーツタイツのCW-Xとか色々と買ってしまった。

お仕事

Courseraで学んだデータ分析・機械学習系の知識を活かして、これまでとちょっと違った方面の仕事に挑戦してみた1年で、川崎グランシティモール実証実験 なんかにも関われて、まだまだこれからだけれど、面白い仕事ができた。

その他にはFMSPLEで発表したりとか。

計算機代数

代数的実数や量化子除去に関して勉強したり簡単な実装をしたりしていたので、思いきって Haskellで計算機代数勉強会という勉強会をゴールデンウィークに開催。結構マニアックなトピックにも関わらず参加者20人以上も集まり、なかなか楽しかった。 その後も、Berlekamp-Zassenhaus アルゴリズムを理解して実装したりとちょっとは進めた。あと、これまでなんとなくしか理解できていなかった sugar flavor とか syzygy とかについて、計算機代数勉強会での@mr_konnさんの発表で理解できて、自分でも実装してみようと思っていたけれど、これは実装しないままになってしまったなぁ……

Max-SAT Evaluation 2013

昨年は Pseudo Boolean Competition 2012 に参加したが、今年は開催されなかったので、代わりに Max-SAT Evaluation 2013 に参加。 今回は自分のtoysatに加えて、GLPKSCIP に簡単なフロントエンドを付けたもので参加。

自作のtoysatの方は、昨年の Pseudo-Boolean Competition 2012 ではそれなりに健闘した(と思ってる)けれど、今回はあまり振るわなかったねぇ。 まあ、色々とやりたいことはあったけれど、結局、目的関数値の線形探索ベースの単純なアルゴリズムのやつで投稿してしまったからなぁ……

SCIPの方は、SCIP本体のバグが見つかって、SCIP開発者に報告・修正してもらったものを再投稿させてもらったりと、意外と大変だったけれど、その甲斐もあって Partial Max-SAT 部門の Crafted ベンチマークで2位に入賞して嬉しかった。(別にフロントエンドを書いただけの自分がすごい訳ではないのだけれど)

英語

英語は今年はあんまり頑張れなかった。 Courseraでは英語を使って学習してはいたし、日本開催の国際会議SPLC2013の併設ワークショップFMSPLE 2013で発表はしたものの、海外には一回も行かなかったし。それに、レアジョブも後半は殆どやらなくなってしまったしなぁ。2014年はその辺りもっと頑張りたいところ。

関連エントリ


2013-12-03 [長年日記]

λ. 自動微分でロジスティック回帰

自動微分を使って線形回帰をしてみる」と「nonlinear-optimizationパッケージで遊ぶ」で線形回帰をして遊んでみたので、次は予定通り自動微分でロジスティック回帰を試してみることに。 例として何を使おうか悩んでいたのだけれど、たまたま ロジスティック回帰で分類を試す - Negative/Positive Thinking (2013-11-29) という記事を(@chezouさんのtweetで)見かけたので、これと同じデータで試してみることにした。

使用するデータ

というわけで、使用するデータはこれ。

コード

実装にはnonlinear-optimizationパッケージを、自動微分を行うadパッケージと組み合わせて便利に使うための、拙作のnonlinear-optimization-adパッケージを使用。

コスト関数としては Courseraの Machine Learning のコースで扱った形の、正則化なしのコスト関数と、L2正則化を行うコスト関数の場合の両方を試す。

import Control.Monad
import Data.Map (Map)
import qualified Data.Map as Map
import Data.Set (Set)
import qualified Data.Set as Set
import qualified Numeric.Optimization.Algorithms.HagerZhang05.AD as CG
import Text.Printf

main :: IO ()
main = do
  trainingData <- liftM parse $ readFile "a9a"
  testData     <- liftM parse $ readFile "a9a.t"

  let features :: Set String
      features = Set.unions [Map.keysSet fs | (_, fs) <- trainingData]

      cost, costReg :: RealFloat a => Map String a -> a
      cost theta    = sum [ if y then - log p else - log (1 - p)
                          | (y,x) <- trainingData, let p = h theta (fmap realToFrac x) ]
                      / fromIntegral (length trainingData)
      costReg theta = cost theta + (lambda / (2 * fromIntegral (length trainingData))) * sum [x**2 | (feat,x) <- Map.toList theta, feat/=""]

      params = CG.defaultParameters
               { CG.printFinal  = True
               , CG.printParams = True
               , CG.verbose     = CG.Verbose
               , CG.maxItersFac = 50 / fromIntegral (1 + Set.size features)
               }
      grad_tol = 0.00001
      theta0 = Map.fromAscList [(feat,0) | feat <- Set.toAscList (Set.insert "" features)]

  (theta1, _result1, _stat1) <- CG.optimize params grad_tol theta0 cost
  (theta2, _result2, _stat2) <- CG.optimize params grad_tol theta0 costReg

  forM_ [("without regularization", theta1), ("with regularization", theta2)] $ \(s, theta) -> do
    printf "[%s]\n" s
    forM_ [("training set", trainingData), ("test set", testData)] $ \(name, samples) -> do
      let correct  = filter (\(y,x) -> predict theta x == y) samples
          accuracy :: Double
          accuracy = fromIntegral (length correct) / fromIntegral (length samples)
      printf "  %s accuracy: %d/%d = %f %%\n" name (length correct) (length samples) (100*accuracy)

-- regularization parameter
lambda :: Num a => a
lambda = 1

-- hypothesis
h :: RealFloat a => Map String a -> (Map String a -> a)
h theta x = sigmoid $ sum [(theta Map.! feat) * val | (feat, val) <- ("", 1) : Map.toList x]

predict :: Map String Double -> (Map String Double -> Bool)
predict theta x = h theta x > 0.5

sigmoid :: RealFloat a => a -> a
sigmoid x = 1 / (1 + exp (-x))

parse :: String -> [(Bool, Map String Double)]
parse = map f . lines
  where
    f l = (readInt w > 0, Map.fromList [(name, read val) | s <- ws, let (name,_:val) = break (':'==) s])
      where
        (w:ws) = words l
        readInt :: String -> Int
        readInt ('+':s) = read s
        readInt s = read s

ちょっと試してみたところなかなか収束しなかったので、参照したブログと同じくイテレーション回数50回で打ち切ることに。これは、nonlinear-optimizationパッケージのoptimize関数は、イテレーション回数の指定がパラメータ数に対する倍率になっていて、ちょっと面倒だった。

結果

結果は以下のとおりで、参照したブログの数字ともほぼ一緒だし、まあこんなものだろう。

[without regularization]
  training set accuracy: 27644/32561 = 84.89911243512176 %
  test set accuracy: 13846/16281 = 85.04391622136232 %
[with regularization]
  training set accuracy: 27648/32561 = 84.91139707011456 %
  test set accuracy: 13850/16281 = 85.06848473680978 %

2013-11-28 [長年日記]

λ. nonlinear-optimizationパッケージで遊ぶ

自動微分を使って線形回帰をしてみる」の続きで、最急降下法を自分で書いてみるかと思ったのだけれど、そういえばHackageにnonlinear-optimizationパッケージというのがあったのを思い出した。以前から、一回試してみたいと思っていたパッケージだったので、せっかくなのでこれを試してみる。 試してみるバージョンは現時点での最新版の0.3.7で、Hager and Zhang *1 のアルゴリズムを同著者らが実装したCG_DESCENTライブラリのバインディングが提供されている。

nonlinear-optimizationパッケージのAPIはadパッケージのgradientDescent関数に比べるとちょっと面倒臭いAPIだけれど、前回のコードを書き換えてみたのが以下のコード。 nonlinear-optimizationパッケージは自動微分とかはしないので、勾配を計算する関数を明示的に渡す必要があり、adパッケージのgrad関数で自動微分した結果の関数を渡すようにしてある。

module Main where

import Control.Monad
import qualified Data.Vector as V
import Numeric.AD
import qualified Numeric.Optimization.Algorithms.HagerZhang05 as CG
import Text.Printf
import qualified Text.CSV as CSV

main :: IO ()
main = do
  Right csv <- CSV.parseCSVFromFile "galton.csv"
  let samples :: [(Double, Double)]
      samples = [(read parent, read child) | [child,parent] <- tail csv]      
      -- hypothesis
      h theta x = (theta V.! 0) + (theta V.! 1) * x
      -- cost function
      cost theta = mse [(realToFrac x, realToFrac y) | (x,y) <- samples] (h theta)
  (theta, result, stat) <-
    CG.optimize
      CG.defaultParameters{ CG.printFinal = True, CG.printParams = True, CG.verbose = CG.Verbose }
      0.0000001                  -- grad_tol
      (V.fromList [0::Double,0]) -- Initial guess
      (CG.VFunction cost)        -- Function to be minimized.
      (CG.VGradient (grad cost)) -- Gradient of the function.
      Nothing                    -- (Optional) Combined function computing both the function and its gradient.
  print theta

-- mean squared error
mse :: Fractional y => [(x,y)] -> (x -> y) -> y
mse samples h = sum [(h x - y)^(2::Int) | (x,y) <- samples] / fromIntegral (length samples)

で、これを実行してみると……

PARAMETERS:

Wolfe line search parameter ..................... delta: 1.000000e-01
Wolfe line search parameter ..................... sigma: 9.000000e-01
decay factor for bracketing interval ............ gamma: 6.600000e-01
growth factor for bracket interval ................ rho: 5.000000e+00
growth factor for bracket interval after nan .. nan_rho: 1.300000e+00
truncation factor for cg beta ..................... eta: 1.000000e-02
perturbation parameter for function value ......... eps: 1.000000e-06
factor for computing average cost .............. Qdecay: 7.000000e-01
relative change in cost to stop quadstep ... QuadCufOff: 1.000000e-12
factor multiplying gradient in stop condition . StopFac: 0.000000e+00
cost change factor, approx Wolfe transition . AWolfeFac: 1.000000e-03
restart cg every restart_fac*n iterations . restart_fac: 1.000000e+00
stop when cost change <= feps*|f| ................. eps: 1.000000e-06
starting guess parameter in first iteration ...... psi0: 1.000000e-02
starting step in first iteration if nonzero ...... step: 0.000000e+00
factor multiply starting guess in quad step ...... psi1: 1.000000e-01
initial guess factor for general iteration ....... psi2: 2.000000e+00
max iterations is n*maxit_fac ............... maxit_fac: 1.797693e+308
max expansions in line search ................. nexpand: 50
max secant iterations in line search .......... nsecant: 50
print level (0 = none, 2 = maximum) ........ PrintLevel: 1
Logical parameters:
    Error estimate for function value is eps
    Use quadratic interpolation step
    Print final cost and statistics
    Print the parameter structure
    Wolfe line search ... switching to approximate Wolfe
    Stopping condition uses initial grad tolerance
    Do not check for decay of cost, debugger is off
iter:     0 f =   4.642373e+03 gnorm =   9.306125e+03 AWolfe =  0
QuadOK:  1 initial a:   1.070618e-04 f0:   4.642373e+03 dphi:  -8.662251e+07
Wolfe line search

iter:     1 f =   5.391563e+00 gnorm =   3.269753e-02 AWolfe =  0
QuadOK:  1 initial a:   1.632104e+01 f0:   5.391563e+00 dphi:  -1.069411e-03
Wolfe line search
RESTART CG

iter:     2 f =   5.387312e+00 gnorm =   1.558689e+01 AWolfe =  1
QuadOK:  0 initial a:   3.264209e+01 f0:   5.387312e+00 dphi:  -2.430188e+02
Approximate Wolfe line search

iter:     3 f =   5.374303e+00 gnorm =   3.196800e-02 AWolfe =  1
QuadOK:  1 initial a:   3.559386e+00 f0:   5.374303e+00 dphi:  -1.022238e-03
Approximate Wolfe line search
RESTART CG

iter:     4 f =   5.288871e+00 gnorm =   2.808049e-02 AWolfe =  1
QuadOK:  1 initial a:   2.376426e+01 f0:   5.288871e+00 dphi:  -7.887345e-04
Approximate Wolfe line search

iter:     5 f =   5.279499e+00 gnorm =   1.301306e+01 AWolfe =  1
QuadOK:  0 initial a:   4.752853e+01 f0:   5.279499e+00 dphi:  -1.693871e+02
Approximate Wolfe line search
RESTART CG

Termination status: 0
Convergence tolerance for gradient satisfied
maximum norm for gradient:  5.513258e-09
function value:             5.000294e+00

cg  iterations:                   6
function evaluations:            15
gradient evaluations:            11
===================================

fromList [23.941530180694265,0.6462905819889333]

今度は一瞬で収束した! やった!

*1  Hager, W. W. and Zhang, H. A new conjugate gradient method with guaranteed descent and an efficient line search. Society of Industrial and Applied Mathematics Journal on Optimization, 16 (2005), 170-192.


2013-11-26 [長年日記]

λ. 自動微分を使って線形回帰をしてみる

Conal Elliott の Beautiful differentiation の論文を読んだ話 の続き。 自動微分の機械学習系での応用は分かるという話で、せっかくなので一番簡単な線形回帰でも使って試してみることに。

データ例

例としては、Rでお馴染みのgaltonのデータを使う。 これは、リンク先にも書いてあるけど、Galtonさんが1885に両親の身長と子供の身長の関係を分析したデータセットで、928サンプルが含まれている。RのUsingRパッケージに含まれているので、これをRで以下のようにして、galton.csvに保存しておく。

> install.packages("UsingR")
> library(UsingR)
> data(galton)
> write.csv(galton, "galton.csv", row.names = FALSE)

せっかくなので、Rで線形回帰をして散布図上にプロットしておく。

> lm <- lm(galton$child ~ galton$parent)
> summary(lm)

Call:
lm(formula = galton$child ~ galton$parent)

Residuals:
    Min      1Q  Median      3Q     Max 
-7.8050 -1.3661  0.0487  1.6339  5.9264 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)   23.94153    2.81088   8.517   <2e-16 ***
galton$parent  0.64629    0.04114  15.711   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 2.239 on 926 degrees of freedom
Multiple R-squared: 0.2105,	Adjusted R-squared: 0.2096 
F-statistic: 246.8 on 1 and 926 DF,  p-value: < 2.2e-16

> plot(galton$parent, galton$child, pch=19, col="blue")
> lines(galton$parent, lm$fitted, col="red", lwd=3)

自動微分のコード

Conal Elliott のコードを使おうかと思ったけれど、Hackageにekmettプロダクトのadパッケージがあったので、それを使うことに。バージョンは3.4を使う。 そういえば、先日開催されていたekmett勉強会でも、@nebutalabという方がadの紹介をしていた。

最急降下法

で、grad関数で自動微分して勾配を求めて、それを使って最急降下法するコードを自分で書こうと思ったら、ちょうど都合の良いことにgradientDescentという関数が用意されていた。 なので、これを使うことにして書いてみたのが以下のコード。 gradientDescentには普通に目的関数costと初期値を引き渡しているだけで、勾配を計算するための関数などを全く引き渡していないことに注目。

-- Galton.hs
module Main where

import Control.Monad
import Numeric.AD
import Text.Printf
import qualified Text.CSV as CSV

main :: IO ()
main = do
  Right csv <- CSV.parseCSVFromFile "galton.csv"
  let samples :: [(Double, Double)]
      samples = [(read parent, read child) | [child,parent] <- tail csv]      
      -- hypothesis
      h [theta0,theta1] x = theta0 + theta1*x
      -- cost function
      cost theta = mse [(realToFrac x, realToFrac y) | (x,y) <- samples] (h theta)
  forM_ (zip [(0::Int)..] (gradientDescent cost [0,0])) $ \(n,theta) -> do
    printf "[%d] cost = %f, theta = %s\n" n (cost theta :: Double) (show theta)

-- mean squared error
mse :: Fractional y => [(x,y)] -> (x -> y) -> y
mse samples h = sum [(h x - y)^(2::Int) | (x,y) <- samples] / fromIntegral (length samples)

実行結果

それで、コンパイルして実行してみると……

[0] cost = 3156.057306315988, theta = [2.6597058526401096e-2,1.8176025390625015]
[1] cost = 2146.158122210389, theta = [4.6848451721258726e-3,0.3193588373256373]
[2] cost = 1459.9671459993929, theta = [2.2758658044906822e-2,1.5543558052476183]
[3] cost = 993.724511817967, theta = [7.872126807824226e-3,0.5363518751316139]
[4] cost = 676.9290401789791, theta = [2.0154694031410444e-2,1.3754888484523264]
[5] cost = 461.6776592383078, theta = [1.0041866689573948e-2,0.6837909525010937]
[6] cost = 315.42191649062323, theta = [1.8389486179011837e-2,1.2539549823963831]
[7] cost = 216.0462832295627, theta = [1.1520222735584755e-2,0.783970560056418]
[8] cost = 148.52403553080242, theta = [1.719418363575862e-2,1.1713769350191798]
[9] cost = 102.64504299700745, theta = [1.252880777988588e-2,0.8520390200495115]
…
[100] cost = 5.391540483678139, theta = [1.5252813280501425e-2,0.9963219150471427]
[101] cost = 5.391540272870362, theta = [1.5259169016306468e-2,0.9963197390947276]
[102] cost = 5.391540062688162, theta = [1.5265580330093717e-2,0.9963213622860531]
[103] cost = 5.3915398529310075, theta = [1.5271945828031475e-2,0.9963198538562846]
[104] cost = 5.391539446704301, theta = [1.5284752349382588e-2,0.996321999765866]
…
[1000] cost = 5.3913214546153885, theta = [2.195100305728909e-2,0.9962239448356897]
[1001] cost = 5.39132107512034, theta = [2.19637098242177e-2,0.9962195158083134]
[1002] cost = 5.391320924102971, theta = [2.1976643022454955e-2,0.996230564935548]
[1003] cost = 5.391320615754077, theta = [2.198280974693417e-2,0.9962155918323486]
[1004] cost = 5.391320339333709, theta = [2.1989373588151114e-2,0.9962277637190161]
…
[10000] cost = 5.389137897180453, theta = [8.88298535971496e-2,0.9952457730583962]
[10001] cost = 5.38913768668423, theta = [8.8836182876215e-2,0.9952431313378857]
[10002] cost = 5.389137477124735, theta = [8.884258017163454e-2,0.995245138984378]
[10003] cost = 5.389137268201876, theta = [8.884892139830691e-2,0.995243314173662]
[10004] cost = 5.389136870256361, theta = [8.886169628000305e-2,0.9952459827149522]

全然収束しない。収束遅すぎだろ。めげずにさらに放置してると……

[100000] cost = 5.367962167905593, theta = [0.7474196350908293,0.9856182795949244]
[100001] cost = 5.367961856713982, theta = [0.7474255920526732,0.9856022065240108]
[100002] cost = 5.367961582358961, theta = [0.7474319755635072,0.9856152902779793]
[100003] cost = 5.367961333033166, theta = [0.747438007468459,0.9856043401613668]
[100004] cost = 5.367961100713857, theta = [0.7474443291979176,0.9856132010817849]
…
[1000000] cost = 5.210315054542271, theta = [6.411550087218479,0.9027479380474824]
[1000001] cost = 5.210314936517216, theta = [6.411554713147153,0.9027442449249086]
[1000002] cost = 5.210314820386824, theta = [6.411559435811994,0.9027471642783936]
[1000003] cost = 5.210314705543812, theta = [6.411564078735267,0.9027446329906037]
[1000004] cost = 5.210314591575562, theta = [6.411568787386962,0.9027465946477329]
…
[5000000] cost = 5.0177281352658305, theta = [18.890815946868038,0.7201815311956332]
[5000001] cost = 5.017728123133143, theta = [18.890817258298373,0.7201790053001258]
[5000002] cost = 5.017728111906337, theta = [18.89081863661461,0.720181051408368]
[5000003] cost = 5.017728101294881, theta = [18.89081995979637,0.7201793288296293]
[5000004] cost = 5.01772809110162, theta = [18.890821328424636,0.7201807127666369]
…
[10000000] cost = 5.001070611475556, theta = [22.875383075131683,0.6618880826205324]
[10000001] cost = 5.00107061095209, theta = [22.87538335250446,0.6618875866914469]
[10000002] cost = 5.001070610463426, theta = [22.875383643001708,0.6618879878893927]
[10000003] cost = 5.001070609998513, theta = [22.875383922680363,0.6618876495886096]
[10000004] cost = 5.001070609549686, theta = [22.87538421127661,0.6618879208540909]

まだかかるのかよ……

[15000000] cost = 5.000328380487849, theta = [23.716478910924902,0.6495830291863548]
[15000001] cost = 5.0003283804652, theta = [23.71647896958202,0.6495829318117019]
[15000002] cost = 5.000328380443848, theta = [23.71647903081446,0.6495830104741416]
[15000003] cost = 5.000328380423461, theta = [23.716479089924043,0.6495829440298111]
[15000004] cost = 5.000328380403642, theta = [23.71647915078346,0.649582997196493]
…
[20000000] cost = 5.0002953079344685, theta = [23.894024471635685,0.6469855785633368]
[20000001] cost = 5.000295307933508, theta = [23.894024484038532,0.64698555944479]
[20000002] cost = 5.000295307932569, theta = [23.89402449694667,0.6469855748657359]
[20000003] cost = 5.0002953079316805, theta = [23.894024509438292,0.6469855618158963]
[20000004] cost = 5.000295307930787, theta = [23.89402452227324,0.6469855722344278]
…
[22987995] cost = 5.000294005873883, theta = [23.922778217862263,0.6465649143875803]
[22987996] cost = 5.000294005873744, theta = [23.922778220364872,0.6465649143543949]
[22987997] cost = 5.000294005873735, theta = [23.922778220366094,0.6465649143543771]
[22987998] cost = 5.000294005873733, theta = [23.922778220367316,0.6465649143543594]
[22987999] cost = 5.000294005873733, theta = [23.922778220367316,0.6465649143543594]
[22988000] cost = 5.000294005873733, theta = [23.922778220367316,0.6465649143543594]

やった! 収束した!

Rだと一瞬だったところ、1.5日くらいかかったし、Rで計算した theta0=23.94153, theta1=0.64629 とはちょっとずれてはいるけれど。

感想

Courseraの Machine Learning のコースで、最急降下法で線形回帰を実装したときは1000ステップくらいで収束してたので、これもそれくらいで収束するかと思ってたんだけれど、こんなにかかるとは以外だった。 gradientDescent関数(ソース)の実装に何か問題があるのか、それともそもそも素朴な最急降下法には荷の重い問題だったのか分からんけれど、2パラメータのこの程度の問題でこれは流石に酷い。

本当はこの後は、ロジスティック回帰とニューラルネットの例も何かやってみようかと思ってたのだけれど、その前にそっちを何とかしないといけないかも。 Machine Learning のコースでも最初は簡単な再急降下法を自前で実装させていたけれど、その後は最適化は予め用意してある実装(fmincg.m)を使うようになってたからなぁ。 まあ、時間かかりすぎて疲れたので、また今度。

【追記】頂いたコメンドなど

【追記】 解決策

本日のツッコミ(全2件) [ツッコミを入れる]

ψ mkotha [入力をおおざっぱに正規化(68を引いてから4で割る)してから最適化したところ、99ステップで(Doubleの最後の桁..]

ψ sakai [おお〜! そういえば、Courseraの講義でも feature scaling が勾配法の収束の速度に効いてくると..]