Home

PuLPで線形計画問題を解く

最近PuLPという線形計画問題 (や混合整数問題) を解くためのPythonライブラリを触ったのでそれの使い方をメモしておきます。

まず線形計画問題について簡単に説明をします。

線形計画問題とは

用語の復習から。

実数$c_1, \dots, c_n \in \mathbb{R}$に対して、次のように与えられる実変数$x_1, \dots, x_n$の関数

線形関数 (linear function) と呼びます。$f$が線形関数で、$b$が実定数のとき、

の関係式を線形等式 (linear equality) といい、

の式を線形不等式 (linear inequality) といいます。

さて線形計画問題 (linear programming problem; LP)とは、有限個の線形等式または線形不等式を満たし、かつ与えられた線形関数を最大化 (最小化) する問題のことです。最大化 (最小化) する線形関数のことを目的関数 (objective function) といい、満たすべき線形等式・線形不等式のことを線形制約式、または線形制約条件 (linear constraint) とよびます。また、線形制約条件すべてを満たす解を実行可能解 (feasible solution) といい、その集合を実行可能領域 (feasible region) といいます。なお、実行可能領域が空集合でない問題を実行可能な問題といい、逆に実行可能領域が空集合である問題を実行不可能な問題であるといいます。1

実行可能解のうち目的関数を最大 (最小) にする解を最適解 (optimal solution) とよびます。

一つ例をあげてみましょう。

例題 (ダイエット問題)

(問題はこちらのサイトのものを数値だけ変えています。)

ある動物園では、飼育しているコアラが太ってしまったため糖質制限をかけた餌を与えてダイエットさせることにしました。コアラの主食であるユーカリには糖質の他にも食物繊維、ミネラル、毒素が含まれていますが、ミネラル以外のものは取りすぎると良くないので、これらにも摂取制限をかけます。 成分構成が異なる2つのユーカリの比を上手く変えることで、なるべく多くのミネラルを摂取させるにはどうすればよいでしょうか?

ユーカリ1gあたりの成分表 (成分の単位: mg)

成分 ユーカリA ユーカリB 摂取上限
糖質 10 14 6400
食物繊維 12 10 6000
毒素 2.5 0.8 900
ミネラル 0.8 0.5

ミネラル摂取量の最大化が目的なので、ユーカリAのグラム数を$x$, ユーカリBのグラム数を$y$とおくと、目的関数である線形関数は、

となります。

また、線形制約条件を式で表すと、

以上の制約条件を満たしながらミネラル摂取量を最大化することが今回の問題のゴールになります。

feasible-region

色がつけられた領域が実行可能領域であり、この中から目的関数を最大化する解を探索することになります。

PuLPによる問題の解法例

PuLPについて

PyPIのページには以下のように書かれています。

PuLP is an LP modeler written in python. PuLP can generate MPS or LP files and call GLPK, COIN CLP/CBC, CPLEX, and GUROBI to solve linear problems.

PuLPはLPを記述するためのモデラーであり、GLPKやCOINなどのソルバーは別に用意してAPIの形で呼び出します。2そのような設計になっているので、最適解を求めるソルバは付け替えができます (デフォルトではCOINというソルバを呼び出す)。

実際にソルバーが最適解を求める際には単体法 (Simplex Method) や内点法 (Internal Point Method) というアルゴリズムがよく使用されます。

それでは、先程の問題をPuLPを使って解いてみましょう。

コード

(ソースコードはgistに置いています。)

PuLP自体はpip install pulpでインストールできます。

まずはユーカリAとユーカリBについて、1mg当たりの糖質、食物繊維、毒素、ミネラルの含有量を変数に代入します。

import pulp

eucalyptus = ['a', 'b']

sugar_content = {
    'a': 10,
    'b': 14
}

fiber_content = {
    'a': 12,
    'b': 10
}

toxin_content = {
    'a': 2.5,
    'b': 0.8
}

mineral_content = {
    'a': 0.8,
    'b': 0.5
}

LpProblem()関数に引数を渡して最適化を解くためのオブジェクトを生成します。 引数の1つ目は問題の名前で、 2つ目は最大化か最小化かを示すためにLpMaximizeLpMinimizeを渡してやります。

prob = pulp.LpProblem('Diet Problem', pulp.LpMaximize)

制約問題の変数を作成します。LpVariable()関数で変数一つずつ作ることもできますが、今回はLpVariable.dicts()で一気に作成をします。lowBoundで変数が取る定義域の下限を指定できます。

eucalyptus_vers = pulp.LpVariable.dicts('Eucalyptus', eucalyptus, lowBound=0)

目的関数を先程作成したprobオブジェクトに追加します。これは(目的関数, 目的関数の説明)のtupleを+演算子でLpProblemのオブジェクトに作用させることで表現できます。今回は0.8x + 0.5yが目的関数なので以下のように表すことができます。

prob += pulp.lpSum([mineral_content[i] * eucalyptus_vers[i] for i in eucalyptus]), 'Total Mineral Content'

制約条件も同じようにprobオブジェクトに+で追加できます。

prob += pulp.lpSum([sugar_content[i] * eucalyptus_vers[i] for i in eucalyptus]) <= 6400, 'Total Sugar Content'
prob += pulp.lpSum([fiber_content[i] * eucalyptus_vers[i] for i in eucalyptus]) <= 6000, 'Total Fiber Content'
prob += pulp.lpSum([toxin_content[i] * eucalyptus_vers[i] for i in eucalyptus]) <= 900, 'Total Toxin Content'

solve()を呼び出すことで最適化をソルバが解いてくれます。

prob.solve()

solve()を呼び出した後、最適解が得られているかは以下のコードを実行すれば確認できます。

pulp.LpStatus[prob.status]

出力が'Optimal'ならば最適解が得られています。

解の出力は以下。

for v in prob.variables():
    print(f'{v.name} = {v.varValue}')

print(prob.objective)
print(f'optimized value: {pulp.value(prob.objective)}')

output:

Eucalyptus_a = 277.03704
Eucalyptus_b = 259.25926
0.8*Eucalyptus_a + 0.5*Eucalyptus_b
optimized value: 351.25926200000004

一応、最適解が得られているか図示して確認しましょう。

def plot_feasible_region(obj_value: float):
    x = np.arange(0,400,0.01)
    y0 = 0 * x
    y1 = - 5 / 7 * x + 1600 / 3
    y2 = - 6 / 5 * x + 600
    y3 = - 25 / 8 * x + 1125
    y_obj = -1.6 * x + 2 * obj_value

    y = np.minimum(np.minimum(y1, y2), y3)
    plt.plot(x, y_obj, color='r', label='目的関数 (ミネラル)')
    plt.plot(x, y1, color='k', label='糖質')
    plt.plot(x, y2, color='b', label='食物繊維')
    plt.plot(x, y3, color='g', label='毒素')
    plt.fill_between(x, y0, y, facecolor='y', alpha=0.5)

    plt.ylim(0, 600)
    plt.legend(prop={'size': 14})

plot_feasible_region(pulp.value(prob.objective))

feasible-region

目的関数が実行可能領域を通る中で$y$切片が一番大きい値になっていることが分かります。たしかに最適な値を出力できていることが確認できました。

補足

問題に等号のない不等式を含む場合は線形計画問題とは呼びません。これは例をみると理由がすぐわかって、以下の簡単な問題

でも、$x_1 = 5$が実行可能領域ではないため最適解を得られない事態が発生するためです。

なお、PuLPは非線形計画問題は解けないので、目的関数や制約条件に非線形な式を入れようとするとエラーを返します。

# 0.2 * x * y == 100 の制約条件
prob += 0.2 * eucalyptus_vers['a'] * eucalyptus_vers['b'] == 100, 'Dummy Non-Linear Constraiant'

出力:

TypeError: Non-constant expressions cannot be multiplied

参考

実行可能な問題でも最適解が存在しない場合もあり、そのような問題は非有界 (unbounded) であるというそうです。それ以外の問題は有界 (bounded) であるといい、LPが最適解を持つことは「問題が実行可能かつ有界」であることと同値のようです。

非線形計画問題を解きたかったらPyomoとか別のパッケージもあるらしい

Read more

平滑化スプライン

スプライン関数

元の現象が連続的なものであっても、実験や統計の結果として得られるデータは離散的なものです。そのため、測定しなかった入力に対する出力が欲しい場合は、手元にある離散的なデータから推定する必要があります。

古典的な推定法としては、Lagrangeの補完法というものがあります。これは、$n$個の点$\{(x_i, y_i)\}_{i=1}^n,\ \ i\neq j \Rightarrow x_i \neq x_j$が与えられたときに多項式を用いてデータを推定するものです。具体的には、

で表される$P(x)$が補間多項式となります。これは全てのデータ点${(x_i, y_i)}_{i=1}^n$を通る、$n-1$次の多項式です。

しかし、このような多項式による近似は良い近似にならない場合があります。例えば以下の例を見てみましょう。

> library('ggplot2')
> library('polynom')
> x <- seq(0, 10, by=1)
> y <- sin(x) + rnorm(11, sd=0.1)
> dat <- data.frame(cbind(x, y))
> poly.calc(x, y)
0.07066721 - 6.503391*x + 18.75596*x^2 - 18.35645*x^3 + 9.239577*x^4 - 2.754634*x^5 + 0.5117172*x^6 - 0.05956099*x^7 +  
0.004201839*x^8 - 0.0001633528*x^9 + 2.667467e-06*x^10 
> f <- function(x) {return(0.07066721 - 6.503391*x + 18.75596*x^2 - 18.35645*x^3 + 9.239577*x^4 - 2.754634*x^5 + 0.5117172*x^6 - 0.05956099*x^7 +  
+ 0.004201839*x^8 - 0.0001633528*x^9 + 2.667467e-06*x^10)}
> ggplot(dat, aes(x=x, y=y)) + geom_point(size=5, col='blue') + stat_function(fun = f, size=1.25, alpha=0.4, xlim=c(-0.3, 10.3))

lagrange_interpolation

上の例はsinカーブにノイズを乗せた点11個からLagrange補間を行った結果を図示したものです。図のように、部分的には元の関数の振舞いを表現できていますが、特にデータの端よりも外にある範囲の補間(補外といいます)が全くおかしなことになっています。これはデータがない部分を表そうというのがそもそも難しいというのもあるのですが、それにしてももう少し上手いこと推定して欲しいところです。

Lagrange補間がデータの端より外にある範囲で変化が激しくなるのは、一つの多項式で全体を表そうとしたことが問題であると考えられます。そこで、単一の多項式によって補完するのではなく、ある範囲のときにだけ多項式になる関数 (区分的多項式) をいくつも用意し、それらを滑らかっぽく繋げることで補完しようというのがスプライン平滑化 (spline smoothing) の気持ちです。

区分多項式とスプライン

ここで用語の定義をしておきます。$a = \xi_0 < \xi_1 < \dots < \xi_n = b$となるような実数の増加列$\{\xi\}$を考え、$[a, b]$上の連続関数で各区間$[\xi_{i-1}, \xi_i]$に制限したときに多項式になるようなものを区分的多項式といいます。区分的多項式$P$は$n$個の多項式$p_0, \dots, p_{n-1}$を用いて、

と定義できます。また、区分的多項式の次数はこのように定義します。

つまり$n$個ある多項式のうちの最大次数を区分的多項式の次数にする、というわけです。

ここで、高々$k$次の区分的多項式で、$(a, b)$の範囲で$C^{k-1}$に属するものを高々$k$次のスプライン関数と呼び、$\{\xi_{i}\}$をノット、$\boldsymbol{\xi} = (\xi_1, \dots, \xi_{n-1})$をノットベクトルと呼びます。 つまり区分的多項式モデルの関数に$k-1$次の微分連続性を求めたものが$k$次のスプラインということです。

さらに、両端の外側に線形制約 (2次以上の項が存在しない) という条件をスプラインに課したものを自然スプラインといいます。$k$次のスプラインに外側に線形制約を課したものは$k$次の自然スプラインですね。

平滑化スプライン

回帰のテクニックとして、RidgeやLassoといった正則化項を導入してモデルが過度に複雑になることを防ぐ方法が知られています。今回はその正則化と同じように、回帰関数の曲率がなるべく小さくなるように関数を求めてみましょう。つまり、

を解いてやればいいことになります。ここで$\lambda$は平滑化パラメータで、$\lambda$が大きいほど曲率が小さい関数が好まれるようになり、$\lambda \to \infty$で$f$は線形モデルになります。

さて、$(1)$の解はどのような関数になるでしょうか?実は、$\{x_i\}$をノットとする3次の自然スプラインになります。これを今から証明していきましょう1。なお、このスプラインを平滑化スプラインと呼びます。

(1)の解が3次自然スプラインになることの証明

$\tilde{g}(x)$を$(1)$の解とする。$\tilde{g}$が各$x_i$において通る点を$z_i = \tilde{g}(x_i)$とおく。一方、$g(x)$を$\{x_i\}$のノットとし、$\{x_i, z_i\}$を通るような3次自然スプラインとする。このような3次自然スプラインは自由度が$n$であることから必ず存在する。

このとき、全ての$x \in [a, b]$において$g(x) = \tilde{g}(x)$であることを示す。

$\tilde{g}(x)$は$(1)$式の解であることから、

を満たす。

次に、$\int_a^b \tilde{g}^{\prime\prime}(x)^2dx \ge \int_a^b g^{\prime\prime}(x)^2dx$を証明する。$h(x) := \tilde{g}(x) - g(x)$という関数$h$を考えると、

$\int_a^b h^{\prime\prime}(x)g^{\prime\prime}(x) dx$について部分積分を行うと、

$g$は自然スプラインなので、$g^{\prime\prime}$は両端では$0$になる。従って$(3)$の第一項は$0$である。

さらに、$g$は3次の区分的多項式なので、$g^{\prime\prime\prime}(x)$は各区間内で定数になる。そこで、区間$(x_j, x_j + 1)$における$g^{\prime\prime\prime}(x)$の値を$c_j$とおくと、

よって$(3)$式の第二項も$0$なので、$\int_a^b h^{\prime\prime}(x)g^{\prime\prime}(x) dx = 0$である。

また、両辺が等しいので$h(x)$が全ての$x \in [a, b]$において$0$のときのみである。$(5)$

$(2)$, $(4)$より$\int_a^b \tilde{g}^{\prime\prime}(x)^2dx = \int_a^b g^{\prime\prime}(x)^2dx$。$(5)$に気をつけると、全ての$x \in [a, b]$について、

(証明終)

なお、$\lambda$は実際にはどうやって求めるかというと、一般化クロスバリデーション (GCV) という手法を使うことで求めるようです (この部分を勉強する余裕ができて理解できれば記事にします。。)。

Rでの実装例

最後にRで平滑化スプラインを行った例を確認してみましょう。今回使用する関数はmgcvというパッケージのgam関数で、gamはGeneralized Additive Model (一般化加法モデル) の略です。一般化加法モデルは平滑化スプラインの説明変数が複数になったバージョンのモデルですが、これを使って平滑化スプラインのモデルも表現できるのでこちらを使用します。

> library('mgcv')
> x <- seq(0, 10, by=1)
> y <- sin(x) + rnorm(11, sd=0.1)
> dat <- data.frame(cbind(x, y))
> gam.model <- gam(y ~ s(x), data=dat)
> plot(gam.model, residuals = T, pch = 19, col='black', xlim=c(-2, 12))

smoothing-spline

図のように、Lagrange補間がデータの端点の外側が11次式の項で急激に変化していたのに対し、平滑化スプラインでの補間は端点の外側が線形であることが保証されているので変化が比較的穏やかになりました。なお、実際のデータ解析では与えられたデータ点全てをノットに使うのではなく、適当な方法で間引いたものをノットとすることが多いようです2

参考文献

1: 統計的学習の基礎 5章前半(~5.6)の証明ほぼそのままです
2: Thin-plate spline法に書いてありました

Read more

Newton Method 2

前回の記事の続き。

関数$f$が狭義凸関数であるときにヘッセ行列が正則になる証明を行っていく。以下では、ヘッセ行列の各成分は実数であるとする。

流れ

先に証明の流れを書いておく。

1. 正定値対称行列が正則であることを証明する

2. ヘッセ行列が正定値対称行列であることを証明する
2.1. 狭義凸関数の性質である

を証明する。

2.2. $f$を2次のテーラー展開して、ヘッセ行列が式中に現れることを確認する。

2.3.下の式が導けるので、 凸関数の性質からヘッセ行列が正定値であることが証明される

正定値対称行列が正則であることの証明

まず正定値対称行列について。

$n\times n$行列$A$が長さ$n$の任意の実数値ベクトル$z$について、$z^TAz > 0$を満たすとき、$A$を正定値行列(positive definite matrix)という1。更に、正定値行列$A$が対称行列であるときその行列$A$は正定値対称行列であるという。

$A$は対称行列であるので、その成分が実数であるとき対角化可能。よって、適当な直交行列$P$と対角行列$D$を用いて、

のように表せる。ここで、任意のベクトル$z (z\in\mathbb{R^n})$を用いると、

のように式変形できる(最後の不等号は正定値行列の定義より)。 $y$が任意のベクトルであることから、$(1)$式より$A$を対角化した行列の対角成分は全て正であることが分かる。

さて、$A$の行列式を考えると、

対角行列の行列式はその対角成分の積であるので、行列式は正。よって正則である。

ヘッセ行列が正定値対称行列であることの証明

狭義凸関数の性質である

を証明する。

狭義凸関数の性質より、

$h:= y - x$とおくと、

これを式変形すると、

が任意の$k\in\mathbb{N}$で成り立つことが分かる。よって微分の定義より、

証明できた。

ところで、この式は感覚的にも分かりやすい式であると思う。

cpp-snip-wrong-result

要するに、狭義凸関数$f$上のある地点を通る接線は$f$以下の値を取る、ということを主張しているだけである。

次に、$f$について2次のテーラー展開を行うと、

上式を満たすような$c$が存在することがわかる。なお、

である。$(3)$式を変形すると、

となる($(2)$式よりこの式は正であることがわかる)。 $(4)$式よりヘッセ行列は正定値行列であることがわかり、よって正則であることが示された。

reference

BASIC PROPERTIES OF CONVEX FUNCTIONS


Definiteness of a matrix - Wikipedia

https://lecture.ecc.u-tokyo.ac.jp/~nkiyono/gocho/en10sol.pdf

Read more

Newton Method

Newton法について勉強したので記録する1。なお、この記事ではベクトルをボールド体にせず、スカラと同じように表記する。

語句の定義

狭義凸関数 (strictly convex function)

凸集合$X$を定義域とする関数$f$が、任意の2点 $x_0, x_1 \in X$, $0 \le \lambda \le 1$である任意の実数$\lambda$に対して、

が成立するとき、$f$は狭義凸関数とよばれる。

(凸集合の定義は前回の記事を参照のこと。)

つまり狭義凸関数は凸関数の定義から等号を除いたバージョンということである。

Newton法がやっていること

Newton法は、狭義凸関数 (strictly convex function) $f$に対して、step wiseに$f$を最小化するパラメータを求めるアルゴリズムである。

なお、以下の説明で出てくる$k$は非負整数であり、$f$を最小化するためにパラメータを更新する回数のことである。

最小化したい$f$は$C^n (n \ge 2)$級の狭義凸関数であるとする。$f(\theta)$について、2次のテーラー展開で近似した$\theta_k$周りの関数を$f_{\text{quad}}(\theta_k)$とすると、 ($\theta \in \mathbb{R}^n, n \in \mathbb{N}$)

ただし、$H_k$は$f_k$についてのヘッセ行列であり、

である。

ここで、以下のように記号を導入する。

$A = \frac{1}{2}H_k$

$b = \nabla f(\theta_k) - H_k\theta_k$

$c = f_k - (\nabla f(\theta_k))^{T}\theta_k + \frac{1}{2}\theta_k^TH_k\theta_k$

このようにおくと、式$(1)$は

のように変形できる。これは$f_{\text{quad}}$を$\theta$について整理した形である。

さて、ここで$f_{\text{quad}}$の極値を求める。

ところで、式$(2)$と、$f$が$C^n (n \ge 2)$級の関数であることから$A$の各成分の微分は順序交換が可能である。よって$A$は対称行列であることがわかるので、

となる。よって、$A$が正則であるならば

上式を満たす$\tilde \theta$は元の関数$f$を二次近似した関数$f_{\text{quad}}$が最小値となるパラメータ (minimizer)なので、$\tilde \theta$は$\theta_k$よりも$f$の最小値に近いことが期待できる。

よって、$k$番目のstep sizeを$\eta_k \in \mathbb{R}$とすると、2

のように$\theta$を更新することを$\theta$の値が収束するまで繰り返せば、$\theta$がminimizerに近い値になっているであろう、というアルゴリズムの考え方である。

ところで、$A$、つまりヘッセ行列が正則である場合はminimizerを求めることができたが、$f$が狭義凸関数であるときにヘッセ行列が正則になる証明はしていない。これは次回以降の記事に回そうと思う (本記事執筆時点で証明はできているはず)。

Reference

  1. Machine Learning: A Probabilistic Perspective, Kevin P. Murphy, 2012
  2. 偏微分の順序交換の十分条件とその証明 高校数学の美しい物語

1: Newton法、 のように値を更新する、という説明が書いてある記事がちらほら見られるけど、これはどのように違うものなのだろうか?違う手法に同じ名前がついてしまっているのか、同じものだけど流派が違うよ、くらいのものなのか。。?

2: なんでstep sizeが出てくるのかよく分かっていない。。そのまま$\tilde \theta$に更新するんじゃダメなんだろうか?

Read more

凸関数の和は凸関数

凸関数の和は凸関数になるが、それの証明を確認したので記録しておく。

以下では簡単化のために集合$X$は有限次元ユークリッド空間$R^m$の部分集合であるとする。

まず用語の定義から。

凸結合

$X$の2つの点$\boldsymbol{x}_0$, $\boldsymbol{x}_1$をとり、$0 \le \lambda \le 1, \lambda \in \mathbb{R}$である$\lambda$を一つ定める。このとき、

で表される点の集合$\boldsymbol{x}$を、$\boldsymbol{x}_0$, $\boldsymbol{x}_1$の凸結合という。

===

この凸結合は2点を結ぶ線分のことです (3点以上に対しても凸結合は定義できる)。

(念の為補足するとこれは、

となることから、$\boldsymbol{x}$は$\boldsymbol{x}_1$を通り$\boldsymbol{x}_1 - \boldsymbol{x}_0$の方向に $|\boldsymbol{x}_1 - \boldsymbol{x}_0|$の$\lambda$倍だけ動く点の集合だから2点を結ぶ線分になる。)

凸集合

集合$X$がそこに含まれる任意の2点から作られる凸結合が$X$に含まれるとき、集合$X$は凸集合であるという。

凸関数

凸集合$X$を定義域とする関数$f$が、任意の2点 $\boldsymbol{x}_0, \boldsymbol{x}_1 \in X$, $0 \le \lambda \le 1$である任意の実数$\lambda$に対して、

が成立するとき、$f$は凸関数とよばれる。

凸関数の和は凸関数

$n$を正の整数として、$f_1, \dots, f_n$が凸関数であるとき、$\sum_{i=1}^{n} f_i$は凸関数になる。

[証明]
$\boldsymbol{x}_0, \boldsymbol{x}_1 \in X$, $0 \le \lambda \le 1$とすると、

Reference

http://web.econ.keio.ac.jp/staff/ito/pdf97/me97fuku.pdf

Read more

macOSで開いた画像ファイルをドラッグ&ドロップする

小ネタだけど地味に便利だったのでメモ。

macOSにおいて、画像ファイルをPreviewなどのアプリで開いた後にその画像をどこかに貼り付けたいときがある。 今まではFinderを開き、画像ファイルが置いてある場所まで移動してファイルをドラッグ&ドロップする、という手順を行っていたが、 Previewで開かれるウィンドウの上にあるファイル名が書いてある隣りにある小さなアイコンはドラッグできるので、 それを直接ドラッグしてペーストしたいところにドロップすればコピペが完了する。

drag-icon

ついでに、ファイル名をクリックするとリネームができる。

rename-file

Previewに限らず、ウィンドウの上にアイコンとファイル名が表示されるタイプのアプリはおそらく全てこの操作ができる。

結構長いことmac使ってたけど、今まで知らなかった。。これは常識なのかな?まあ自分のように知らない人にはよい情報だと思うので共有。

Read more

MacにLightGBMを再インストール

LightGBMを使おうと思いpyhtonのnotebook上でimportしようとしたら、エラーが出た。

from LightGBM import LGBMClassifier

Output (Excerpt):

OSError: dlopen(/usr/local/lib/python3.6/site-packages/LightGBM/lib_LightGBM.so, 6): Library not loaded: /usr/local/opt/gcc/lib/gcc/7/libgomp.1.dylib
  Referenced from: /usr/local/lib/python3.6/site-packages/LightGBM/lib_LightGBM.so
  Reason: image not found

この前までは正常に使えていたと思ったのだが。エラーを見るとgccの7系のライブラリを見に行っているようだが、自分の環境には7系はなかった。じゃあなんで今まで動いていたんだろう?8にアップデートしたのが最近だったんだろうか?

と思ってbrewのlogを見てみたら、5/7にgccのバージョンを7.3.0_1 -> 8.1.0に上げていた。なるほど。

エラーの原因もわかったので、動くようにする。一旦LightGBMをアンインストールして再度インストールすることにした。

LightGBMのインストールガイドを見ると、LightGBMはコンパイルにOpenMPを使っており、これはApple Clangじゃサポートしてないからgcc(の7系)を使ってコンパイルしてくれと言っている。gccは8系なので、7系をインストールする。

$ brew install gcc@7

環境:

$ sw_vers
ProductName:  Mac OS X
ProductVersion: 10.13.4
BuildVersion: 17E202

インストールガイドに沿って、cloneしてビルド、インストールする。

$ cd /tmp
$ git clone --recursive https://github.com/Microsoft/LightGBM
$ cd LightGBM
$ export CXX=g++-7 CC=gcc-7
$ mkdir build
$ cd build
$ cmake ..
$ make -j4

LightGBMのインストールはできたが、元々のエラーは/usr/local/opt/gcc/lib/gcc/7/libgomp.1.dylibがないということだったので、パスを確認してみる。

$ ls /usr/local/opt
lrwxr-xr-x 19 fhiyo  7 May 12:00 gcc -> ../Cellar/gcc/8.1.0
lrwxr-xr-x 23 fhiyo  4 Apr 15:57 gcc@5 -> ../Cellar/gcc@5/5.5.0_2
lrwxr-xr-x 21 fhiyo 19 May 18:08 gcc@7 -> ../Cellar/gcc@7/7.3.0

とりあえずこれ書き換えるしかないっぽいのでそうする。

$ mv gcc gcc@8
$ ln -s gcc@7 gcc

これで動くようにはなったけど嫌だなぁ

ちなみに

Issueが4日前に上がっていた。LightGBM and gcc 8 in MacOS: Library not loaded: /usr/local/opt/gcc/lib/gcc/7/libgomp.1.dylib · Issue #1369 · Microsoft/LightGBM

PyPIの現時点の最新である2.1.1はgcc-7になっているようだが、8を使うように変更したプルリクは既に上がっているとのこと。問題が解決するまでそんなには時間かからなさそう?

追記 (2018/05/21)

$ R
dyld: Library not loaded: /usr/local/opt/gcc/lib/gcc/8/libgfortran.5.dylib
  Referenced from: /usr/local/Cellar/r/3.5.0_1/lib/libR.dylib
  Reason: image not found
zsh: abort      R

やっぱり他に影響があった。

Read more

Juliaを使ってKaggleのhouse pricesの問題を試してみた

今回はJuliaという比較的新しめの言語でKaggleをやってみることにした. まだKaggleも大してやれていないのになぜ今まで触ったこと無い言語を試してみたかというと,社内のハッカソンでやることになったから.

Juliaは実行時コンパイルすることでPythonより高速に計算できることを売りにしている言語のようだ (だが今回の簡単なKernelではその良さは享受できず...).どんな言語かはこちらとかを見ればわかると思う→Python使いをJuliaに引き込むサンプル集

題材にしたKaggleのコンペはHouse Prices: Advanced Regression Techniquesで,アイオワ州のエイムズにある家についての様々な情報 (特徴量の数が79個) から家の値段を予測する,というもの.要するに家の値段について回帰をしろ,という問題.

試してみた結果はこちらのGistにKernelを公開している.本当はKaggle上で公開したかったのだが,現時点ではJuliaのKernel公開用の環境が整っていないみたい (昔のversionのKernelは投稿できたようだ).

Juliaのライブラリはまだあまり成熟していないのか,versionの依存関係によって動かなくなる関数があったりwarningが出たりすることが多く (dependency hellというやつ?),中々苦戦した. (まだ見れていないが,Pkg3というパッケージ管理システムを使えば問題が解決されるのかも?) またJuliaは動的型付けながらも,コンパイルされる関係からか型にうるさい印象だった. Pythonではよしなにやってくれていた処理が全然通らず,NAという欠損値を表す特別な値が型チェックですぐに引っかかったり,GridSearchCVがなぜか全然動かなかったりでハッカソンの一日では前処理をまともに行うことができずじまいだった.

言語自体にはLispのようにマクロが使えたり,関数型言語としての側面もあるので慣れれば楽しそうではあるものの,データ分析でぱっと使うにはライブラリが整っていないためまだまだ厳しいという印象.ネット上のサンプルが0.4とか昔のversionのものが多くてコピペで動作確認ができないし,helpで呼び出せるリファレンスもまだ充実していない感じなので勉強する環境が十分に整っていないのでは,と思う.まあ一日ちょっと触ってみただけなので,もう少し触っていればやりやすくなるのかもしれないが.

KaggleのKernelを見てやり方を勉強する方法はとてもよいと思うので,ぜひJuliaのKernelも早く公開できるようになってほしい.今の状態でコード書くのは自分で試行錯誤しなければならない部分が多くて無駄に時間を食ってしまった.

Read more