画像認識について知りたいなら『画像認識』(著:原田達也先生)がお勧め!

画像認識について勉強したかったため、機械学習プロフェショナルシリーズの『画像認識』を読んでみました。300ページ近くありますので、他の機械学習プロフェショナルシリーズより、少し長めです。しかしその分、画像認識の様々な手法が網羅されている印象です。読み終わった後も、辞書的に置いておくような使い方もできそうですね。

読むのにかかった時間はだいたい20時間ぐらい。ただ、理解を深めるために合わせて読んだ他の本や記事、観た動画などの時間も含めるとさらに10時間ぐらいでしょうか。恐らく読者の事前知識によって大幅に変動しそうですが、画像について勉強したことのある人であれば、さらに速く読みこなせるかと思います。

私は画像について本格的に勉強するのは今回が初めてでしたので、他の画像本、ビジョン本との比較してどうなのか、という観点からはコメントできませんが、せっかくですので、内容とその感想をご紹介していきます。

本のリンクはこちら

第1章 画像認識の概要

内容は画像認識の概要から始まります。こちらは画像認識の歴史を概観し、「クラス認識」「物体認識」「インスタンス認識」など様々なタスクが紹介されていきます。中でも「クラス認識」の流れが詳細に解説されます。その流れとは、

  • 入力画像が与えられる
  • サンプリングと記述を行う
  • 統計的特徴抽出を行う
  • コーディング
  • プーリング
  • 分類器を学習する
  • クラスを予測する

となります。これらを一つ一つ、第2章から第5章まで学んでいきます。そしてこのパイプラインを多段に重ねて、最後に分類器を組み入れるEnd to Endなシステムを深い構造と呼び、その話を第6章で学びます。

第2章 局所特徴

画像認識のはじめの処理である局所特徴を学びます。画像の局所的な小領域に着目して、その内容を記述する必要があります。前半が「検出」で、後半が「記述」になりますが、それぞれについて様々な手法があります。フィルターとは何なのか、なぜ画像を微分するのか、という基本から学び、エッジ検出器、コーナ検出器、ブロブ検出器を見ていきます。記述子としてもSIFT、HOG、SURF、BRIEF、BRISK、GISTなどの代表的なものを見ていきます。

局所特徴を勉強すること自体が初めてでしたので、「セルに分けて勾配のヒストグラムを作って、…」のあたりで少し時間がかかりました。Youtubeに上がっているビジョン系の動画を合わせて勉強すると、ビジュアルでスッと頭に入ってきましたので、動画系を合わせることもお勧めです。著者の先生の公開スライドも理解の助けになりました。

第3章 統計的特徴抽出

「実際の画像には、なんらかの外乱やノイズが加わる」ため、「外乱などの影響を受けている画像から得られた局所特徴をそのまま利用」するわけにはいきません。そこで統計的な構造にもとづいて特徴を抽出する話を学びます。主成分分析、フィッシャー線形判別分析、正準相関分析、偏最小2乗法が出てきます。統計学機械学習で学習済みの方であれば、すぐに終わります。

第4章 コーディングとプーリング

「局所特徴を認識に有効な次元数のベクトルに変換する操作を、コーディング」と呼び、「画像領域内に存在する複数のコーディング後の特徴ベクトルを1本のベクトルにまとめる操作を、プーリング」と呼びます。学ぶ手法としては、

  • Bag of Visual Words
  • コードブック
  • GMMスーパーベクトル
  • フィッシャーベクトル
  • VLAD
  • 局所座標符号化
  • ランダムフーリエ盗聴
  • 空間ピラミッド

などです。他にも色々出てきます。

第5章 分類

主に教師あり分類の話です。

分類とは何か、ベイズ決定則、識別関数、最適化の話題から始まり、後半では各手法(パーセプトロン、Adaline、SVM、ロジスティック回帰、ソフトマックス回帰、局所学習、集団学習)が解説されます。最後に分類結果の評価指標(precisionやrecall、accuracy)、交差検証法が出てきます。

第6章 畳み込みニューラルネットワーク

パーセプトロンから始まり、多層ニューラルネットワークへと拡張していきます。この章で必ず理解しておきたいことは、逆伝播と畳込みのことろでしょう。私は同じくMLPシリーズの『深層学習』を読んで理解を深めました。手を動かすのが億劫になったら、Tesla AI DirectorのAndrej Karpathyのブログ記事などでモチベーションを上げたりしました。こちらのブログで紹介されているStanfordのCS231nの動画も、かなり直感的に逆伝播が解説されており、助かりました。

それ以外にディープを理解する上で必須の内容(活性化関数、最適化方法)も学びます。AlexNet、VGGNet、GoogLeNet、ResNetも解説あり。

第7章 物体検出

画像中から人や車などが写っている領域を特定し、そこを四角い領域で囲む物体検出について学びます。最もナイーブな種法であるスライディングウィンドウ法から始まり、様々なテクニック・手法を学んでいきます。後半ではCNNを利用した物体検出としてR-CNN(CVPR 2014)、Fast R-CNN(ICCV 2015)、Faster R-CNN(NIPS 2015)を学びます。このあたりの物体検出について理解することが今回この本を勉強した理由の一つでした。面白い!!

第8章 インスタンス認識と検索

データベース内の画像をすばやく探す画像検索問題。スペクトラルハッシング、K-meansハッシングなど色々と学びます。

第9章 さらなる話題

セマンティックセグメンテーション、画像からのキャプション生成、画像生成と敵対的生成ネットワーク(GAN)などが出てきます。これからも色んな話題が出てきそうですね。

終わりに

ということで、それぞれの章の内容を簡単にまとめてみました。著者の先生の公開スライドも役立ちました。こちらです:

DeepMindのSilver氏の強化学習の講義動画シリーズを観てみた

強化学習の勉強を始めるために、DeepMindのDavid Silver氏の強化学習の講義を録画した動画を観てみました。「強化学習の勉強をするならコレだ!」と色んなところでお勧めされていたので、ずっと気になっていました。

UCLの強化学習の講義を録画したもので、全部で10回。Youtubeに上がっているので無料で勉強することができます。各講義が1時間30分〜1時間40分ぐらいですので、全体の負担もそんなに多くはありません。(通常、大学の講義動画系はこれの倍以上の時間がかかる傾向にある気がします。)

www.youtube.com

強化学習の基礎の基礎から始めて、DQNなど(少なくとも講義録画時点では)ナウい話も出てきます。求められる前提知識は数学(確率、統計、線形代数微積分)と、機械学習(特に教師あり学習)の知識あたりでしょうか。私は強化学習の知識がほとんどゼロの状態から始めて、時々、本なども読みながら動画を進めていきました。Silver氏の授業が非常にわかりやすいので、わからないところがあっても2度、3度聞けばかなり頭に入ってきます。

指定教科書はこちら:

http://ufal.mff.cuni.cz/~straka/courses/npfl114/2016/sutton-bookdraft2016sep.pdf

こちらのリンクは現在準備中の新しいエディションのものです。概ねこの本の流れに沿って授業が進んでいきますので、わからない箇所があればこちらを参照するのもお勧めです。

10回の授業はそれぞれ次の通りです。

Lecture 1: Introduction to Reinforcement Learning
Lecture 2: Markov Decision Processes
Lecture 3: Planning by Dynamic Programming
Lecture 4: Model-Free Prediction
Lecture 5: Model-Free Control
Lecture 6: Value Function Approximation
Lecture 7: Policy Gradient Methods
Lecture 8: Integrating Learning and Planning
Lecture 9: Exploration and Exploitation
Lecture 10: Case Study: RL in Classic Games

第1回目はお話的なイントロ、最終回はゲームへの応用のお話。特に重要なのは第2回〜第7回という印象です。わからないところがあれば、理解できるまで何度も巻き戻して観ることをお勧めします。あとで同じ手法や概念が何度も登場するからです。

授業の講義ページはこちら。授業スライドなどもこちらからダウンロードできます。
http://www0.cs.ucl.ac.uk/staff/d.silver/web/Teaching.html

Youtubeだと再生回数も見れるのが面白いですよね。第1回目はなんと20万人近く観ていることがわかります。ですが第2回目となると8万人、第3回目は5万人、そして最終回まで頑張った人は1万人近くまで激減(笑)。最後まで粘ると満足感を得られますね!

Chainerにサクッと親しむ『Chainerによる実践深層学習』がオススメ

Chainerによる実践深層学習を読んでみました。Chainerにサクッと親しむためには非常にコスパの良い本でしたので、ここで少し内容をまとめておきたいと思います。

内容としては、Chainerの使い方(NumPyの復習から始まり、オブジェクトやChainクラス、Irisを用いた具体例)を紹介した上で、具体的な深層学習手法を実装していく流れになっています。紹介されている手法はオートエンコーダ、word2vec、Recurrent Neural Network及びLSTM、翻訳モデルです。終盤ではCaffeのモデルをChainerから利用する方法も紹介されておりました。そして最終章はGPUの利用ということで、著者の今までのGPUにまつわる経験やCUDA、CuDNN、CuPyなどの話が一通り乗っていてます。

求められる前提知識は、Pythonをデータ分析で使えること、『深層学習(機械学習プロフェッショナルシリーズ)』程度の知識があること、と感じました。載っているコードはすべて公式サイトからダウンロードすることができ、しかもsklearnのデータを用いることが多いのが利点です。なぜなら、データ分析の本を買っても、前処理で躓いてなかなか前に進めないということがたまにあるからです。(私だけでしょうか…。)今回のサンプルコードは、発売されたばかりということでバージョンの差分がゼロor小さいということもあるかもしれませんが、どれも初期状態で実行するだけでちゃんと動きました。

他の深層学習プログラミング本との違いは、解説するフレームワークをChainer一本に絞り、さらに解説する手法も絞っていることかもしれません。そのため、なかなか取り組みやすいです。私は土日に各5時間ぐらいでさっと一通り読んでみましたが、Chainerのイメージが湧いたので、読んでみてよかったですね。(ただし、パラメータの学習には数十時間かかる箇所もありましたので、その時間は除いています。)

交差検証法の説明と実装

良いモデルを選択したい!

与えられたデータに対して学習アルゴリズムを適用しようとするときに、どのようなモデルを選択するべきなのか、ということを考えなければなりません。例えば線形モデルを適用する場合は基底関数の種類や数を決めないといけないですし、過適合を防ぐために正則化項をモデルに付加している場合は正則化パラメータも調整しなければなりません。それでは、どのようにパラメータを調整したり、モデルの種類を選択すればよいのでしょうか?

そもそも私達の目標は、現在手元にあるデータセットを上手く説明できるモデルを見つけることではなく(そういうケースもあるかもしれませんが、一般的には)将来の未知なデータを予測することでしょう。これを汎化と呼び、将来データに対する予測誤差を汎化誤差と呼ぶことが出来ます。これは訓練標本に対する誤差とは異なることに注意が必要です。現在与えられている訓練標本に対する誤差を小さくするためには、パラメータの値を色々と変更し極端な値にしたり複雑なモデルを選択することでいくらでも小さく出来てしまうケースがあります。これは過適合を起こしてしまいますので好ましくありませんね。

交差検証法

とは言っても、将来データがもし現時点で与えられていれば、それをテスト標本とし、汎化誤差を求めることができますが、現時点ではそれはもちろん与えられていません(もしくはモデルを選択する上でテスト標本を選んではいけない設定を考えることが多いです)。そこで、学習データを幾つかの組に分割し、一部をモデル選択・パラメータ調整のためのデータセット(学習データ)とし、残りを汎化誤差を測るためだけに利用します(評価データ)。最もわかりやすいのは、与えられたデータを2つに分割し、片方で学習し、もう片方で誤差を評価します。これをホールドアウト検証と言います。

ホールドアウト検証では、学習データと評価データに分けてしまい、学習データの数が減ってしまいました。また、分割した2つのデータセットのうち、どちらを学習データとし、どちらを評価データとするかによって、結果が大きく変わってしまうこともあります。そこで、さらに洗練された方法として交差検証法があります。

交差検証法では、与えられているデータセットを2個以上のグループに別けます。例えば50サンプルを10サンプル毎の5グループに分けるような例が考えられます。

  • グループ1を評価用として、その他のグループで学習する
  • グループ2を評価用として、その他のグループで学習する
  • グループ3を評価用として、その他のグループで学習する
  • グループ4を評価用として、その他のグループで学習する
  • グループ5を評価用として、その他のグループで学習する

これによって、5つの誤差が計算でき、それらの平均値を今回のモデルの汎化誤差とします。このような場合を5-fold cross validationと呼びます。実際には10-foldが頻繁に用いられます。はじめの話に戻ると、それぞれのモデルに対して交差検証を行い、最も汎化誤差が小さいモデルを採用すべき、という考え方を取ります。

実装しよう!

それでは以下では具体例として、50個のサンプル標本を基に、交差検証法を行います。ここでは、ガウスカーネルモデルと{l2}ノルム正則化のモデルを利用し、ガウスカーネルのバンド幅{h}正則化パラメータ{\lambda}をどのような値に調整すればよいかを調べていきましょう。分割数は5と、4グループ(40サンプル)を用いて学習し、残りの1グループ(10サンプル)で評価します。これをすべてのグループに繰り返し行います。分割方法はランダムに行う必要があります。詳しくは実装コードを御覧ください。

各実験では、パラメータの値を変更した回数分だけ行います。パラメータを2つ調整するため、少し実装が複雑になります。それぞれのパラメータは{0.0001,0.001,0.01,0.1,1,10,100,1000,10000}の中から選択するようにしています。ですので実験数は5×8×8になり、実装ではそれらが3つの"for"に対応しています。

交差検証用のパッケージもありますが、ここでは理解を深めるために一歩一歩実装していきます。

# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt
import random

s = 50
N = 1000
x = np.matrix(np.linspace(-3,3,s)).T
X = np.matrix(np.linspace(-3,3,N)).T
pix = np.pi * x
y = np.cos(np.pi*x)+x+np.random.randn(s,1)

plt.plot(x,y,'bo',label='data samples')
plt.plot(X,np.cos(np.pi*X)+X,'r')

#preparing for cross validation
hhs = 2*np.power([0.0001,0.001,0.01,0.1,1,10,100,1000,10000],2)
ls = np.array([0.0001,0.001,0.01,0.1,1,10,100,1000,10000])
m = 5
randomlist = range(s)
random.shuffle(randomlist)
randomlist = np.array(randomlist) % 5
g=np.zeros([5,len(hhs),len(ls)])
xx = np.tile(np.power(x,2).T,(s,1))+np.tile(np.power(x,2),(1,s))-2*((x).dot(x.T))

#Cross Validation
for i in range(len(hhs)):
    K = np.exp(-xx/hhs[i])      
    for q in range(m):
        ki=K[randomlist!=q,:]
        kc=K[randomlist==q,:]
        yi=y[randomlist!=q].T
        yc=y[randomlist==q].T
        for j in range(len(ls)):
            n = np.sum(randomlist!=q)
            x2 = np.power(x[randomlist!=q],2)
            t = (np.linalg.inv(ki.T.dot(ki)+ls[j]*np.eye(s))).dot(ki.T.dot(yi.T))
            f = kc.dot(t)
            g[q,i,j] = np.mean(np.power((f-yc.T),2))

G = np.mean(g,0)
smallest = np.unravel_index(G.argmin(), G.shape)
h = hhs[smallest[0]]
l = ls[smallest[1]]

交差検証法で選択されたパラメータはバンド幅が{2\times 1^2}正則化パラメータは0.01となりました。それらの値を用いて、関数を図示していきましょう。

XX = np.tile(np.power(x,2).T,(N,1))+np.tile(np.power(X,2),(1,s))-2*((X).dot(x.T))
kk = np.exp(-xx/h)
KK = np.exp(-XX/h)
TT = (np.linalg.inv(kk.T.dot(kk)+l*np.eye(s))).dot(kk.T.dot(y))
F = KK.dot(TT)

#Plot the graphs
plt.cla()
plt.plot(x,y,'bo',label='data samples')
plt.plot(X,np.cos(np.pi*X)+X,'r',label='y=cos(pi*x)+x')
plt.plot(X,F,'g-',label='After cross validation')
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)

下記が交差検証法で得られたパラメータを用いたプロットです。赤線と緑線が近しいことが視覚的に確認でき、ある程度成功しているのではないかと考えることが出来ますね。
f:id:decompose:20160528004725p:plain

さらに勉強するには

今回の記事を作成する上では以下の書籍が非常に参考になりました。是非覗いてみてください。
『イラストで学ぶ機械学習:最小二乗法による識別モデル学習を中心に』(第4章)
『はじめてのパターン認識』(第2章)
『Introduction to Statistical Learning』(第5章)

K平均法によるクラスタリング

 今回はK平均法を紹介します。
 『はじめてのパターン認識』の表記を参考に話を進めていくと、{d}次元の{N}個のデータ{\mathcal{D}=\{x_1,\ldots,x_N\}}があり、{x_i\in\mathbb{R}^d}とします。この時、データの類似度を基に、{K}個のクラスタに分類する方法です。{K}は予め決まったクラスタ数で、データの特性から判断しモデルに与えます。そして各クラスタの代表ベクトルの集合を{\mathcal{M}=\{\mu_1,\ldots,\mu_K\}}とし、{i}番目のサンプルが{k}番目の代表ベクトルが支配するクラスタに帰属するか否かを帰属変数{q_{ik}}(帰属すれば1、そうでなければ0)と表すとすると、K平均法の評価関数を次のように定式化できます:
{J(q_{ik},\mu_k)=\sum^N_{i=1}\sum^K_{k=1} q_{ik}\|x_i-\mu_k\|^2}
これを{\mu_k}に関して最適化してみましょう。
{\frac{\partial J(q_{ik},\mu_k)}{\partial \mu_k}}=2\sum^N_{i=1}q_{ik}(x_i-\mu_k)=0
よって、
{u_k=\frac{\sum^N_{i=1}q_{ik}x_i}{\sum^N_{i=1}q_{ik}}}
クラスタの代表ベクトルは、帰属するデータの平均ベクトルになることがわかります。ただし、{q_{ik},\mu_k}を同時に最適化するのは難しいので、K平均法アルゴリズムによって逐次最適化を行うことが一般的です。

クラスタ数が3、各クラスのデータ数を100、データの次元を2とした場合について、アルゴリズムPythonにて実装すると、次のようになります。

# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt

mean1 = [-5,-5]
mean2 = [0,7]
mean3 = [5,-5]
cov = [[3.5,0],[0,3.5]]

instances = 100

#generate 2-dimensional data for three classes
x = np.random.multivariate_normal(mean1,cov,instances)
y = np.random.multivariate_normal(mean2,cov,instances)
z = np.random.multivariate_normal(mean3,cov,instances)
X = np.vstack((x,y))
X = np.vstack((X,z))

closest = np.floor(np.random.rand(instances*3)*3)
iteration = 10000

#clustering plot in the first iteration
f, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, sharex='col', sharey='row')
ax1.plot(X[closest==0,0],X[closest==0,1],'bo')
ax1.plot(X[closest==1,0],X[closest==1,1],'ro')
ax1.plot(X[closest==2,0],X[closest==2,1],'go')
ax1.set_title('1st iteration')

for i in range(iteration):
    mu1 = np.mean(X[closest==0,:],axis=0)
    mu2 = np.mean(X[closest==1,:],axis=0)
    mu3 = np.mean(X[closest==2,:],axis=0)
    dist = np.zeros((instances*3,3))
    for j in range(instances*3):
        dist[j,0] = np.linalg.norm(X[j,:] - mu1)
        dist[j,1] = np.linalg.norm(X[j,:] - mu2)
        dist[j,2] = np.linalg.norm(X[j,:] - mu3)
    if np.linalg.norm(closest - np.argmin(dist,axis=1)) == 0:
        break
    closest = np.argmin(dist,axis=1)
    if i ==1:
        ax2.plot(X[closest==0,0],X[closest==0,1],'bo')
        ax2.plot(X[closest==1,0],X[closest==1,1],'ro')
        ax2.plot(X[closest==2,0],X[closest==2,1],'go')
        ax2.set_title('2nd iteration')
    if i ==2:
        ax3.plot(X[closest==0,0],X[closest==0,1],'bo')
        ax3.plot(X[closest==1,0],X[closest==1,1],'ro')
        ax3.plot(X[closest==2,0],X[closest==2,1],'go')
        ax3.set_title('3rd iteration')        

ax4.plot(x[:,0],x[:,1],'bo')
ax4.plot(y[:,0],y[:,1],'ro')
ax4.plot(z[:,0],z[:,1],'go')
ax4.set_title('After convergence')

図示をすると、以下のようになります。
f:id:decompose:20160526001019p:plain

Pythonで0〜9の手書き数字画像データに慣れる!

今回は、手書き数字の画像データを扱うことに慣れることを目標にした記事です。データの読み込みや表示、各ラベルの平均値の取得・表示、PCAの実施と可視化などを取り上げました。まずは、sklearn.datasetsからload_digitsをインポートして、図を作成しましょう。

#coding: utf-8
import numpy as np
import pylab
from sklearn.datasets import load_digits
from sklearn.decomposition import PCA
import matplotlib.pyplot as plt

#load digits data
digits = load_digits()

#transform each 8*8 images into a vector of length 64
X = digits.images.reshape((digits.images.shape[0],-1))
Y = digits.target
target_names = digits.target_names

# Show digits
for index, (image, label) in enumerate(zip(digits.images, digits.target)[:10]):
	pylab.subplot(2, 5, index+1)
	pylab.axis('off')
	pylab.imshow(image, cmap=pylab.cm.gray_r, interpolation='none')
	pylab.title('%i' % label)
pylab.show()

f:id:decompose:20160507003300p:plain
次に、各ラベルの画像データの平均値を取得してみます。こちらの場合でも最後に可視化を行います。

meanmatrix = np.zeros((10,8,8))
for i in range(10):
    test = X[y==i].mean(axis=0)
    meanmatrix[i,:,:] = np.reshape(test,(8,8))
    pylab.subplot(2,5,i+1)
    pylab.axis('off')
    pylab.imshow(meanmatrix[i,:,:],cmap=pylab.cm.gray_r, interpolation='none')
pylab.show()

f:id:decompose:20160507222844p:plain
標本数は2000近くありますが、それらを可視化してみます。そのためのPCAで2次元まで落としてみます。

#PCA on load_digits
pca = PCA(n_components=2)
X_r=pca.fit(X).transform(X)

#Visualize
colors = [plt.cm.nipy_spectral(i/10.,1) for i in range(10)]
plt.figure()
for c, target_name in zip(colors, target_names):
    plt.scatter(X_r[y==target_name,0],X_r[y==target_name,1],c=c,label=target_name)
plt.legend(scatterpoints=1)
plt.title('PCA of load_digits')
plt.show()

f:id:decompose:20160507003306p:plain

参考にしたサイト:
Scikit-learnでPCA - Qiita
多層パーセプトロンで手書き数字認識 - 人工知能に関する断創録

Pythonで勾配法を実装してみました。

今回はこの本に出てくる勾配法を実装します。
www.amazon.co.jp

1変数の場合

第3章のアルゴリズム3.1(1変数の勾配法)を、Pythonで実装してみました。gradient_descentは比較のために書いています。本に登場するのはgradient_descent_2の方です。

import numpy as np

def f(x):
    return - x ** 2

def f_prime(x):
    return  - 2 * x

def gradient_descent(x, step, prime, iteration, e):
    for i in range(iteration):
        x_new = x + prime(x) * step
        if abs(f(x_new) - f(x)) < e:
            break
        x = x_new
    return(x,i)

def gradient_descent_2(x, h, prime, iteration, e):
    for i in range(iteration):
        h = np.sign(prime(x)) * abs(h)
        X = x
        X_new = x + h
        
        if f(X) < f(X_new):
            while f(X) < f(X_new):
                if f(X) >= f(X + h):
                    break
                h = 2 * h
                X = X_new
                X_new = X + h
            x = X
            h = h * 0.5
        else:
            while f(X) > f(X_new):
                if f(X) <= f(X - h):
                    break
                h = h * 0.5
                X_new = X_new - h
            x = X_new
            h = 2 * h
        
        if abs(prime(x)) <= e:
            break
            
    return(x,i)

print( gradient_descent(-100,0.1,f_prime,1000,0.0001) )
print( gradient_descent_2(-100,0.1,f_prime,1000,0.0001) )

多変数の場合

次に、多変数の場合で実装しました。下記は2変数の場合に対応しています。この記事(Rで勾配法(勾配降下法・最急降下法)を実装 - Qiita)をかなり参考にさせて頂きました。表記なども基本的には合わせています。難易度はやはり1変数の場合と比べると難しいですね。教科書ではアルゴリズム3.2(多変数の勾配法)です。

import numpy as np

def f(x):
	return -x[0]**2 - 4 * x[1] **2 + 5 * x[0]

def dx(x):
	return -2*x[0] + 5

def dy(x):
	return -8 * x[1]

def search(x0,delta,f,df):
	x = x0
	h = 0.5
	iteration = 10000
	eps = 0.000001
	for i in range(iteration):
		oldx = x
		h = np.sign(df(x)) * abs(h)
		newx = x + h * delta
		if f(x) < f(newx):
			while f(x) < f(newx):
				if f(newx) >= f(newx+h*delta):
					break
				h = 2 * h
				x = newx
				newx = x + h * delta
			h = h * 0.5
		else:
			while f(x) > f(newx):
				if f(x) <= f(x + h * delta):
					break
				h = h * 0.5
				newx = x + h * delta
			x = newx
			h = 2 * h
		if abs(h) < eps:
			break
	return x

def hillclimb(x,f,dx,dy):
	eps = 0.00000001
	iteration = 10000
	for i in range(iteration):
		gradientf0 = np.array([dx(x),dy(x)])
		print('hillclimb i=', i, 'f(x)=', f(x), 'gradientf0=', gradientf0)
		def df(x):
			return np.dot( np.array([dx(x),dy(x)])  ,  gradientf0 )
		print df(x)
		delta = ( gradientf0 ) / np.linalg.norm(gradientf0)
		newx = search(x,delta,f,df)
		diff = newx -x
		x = newx
		if sum(diff * diff) < eps:
			break
		print(delta, sum(diff*diff))
	return (x)


print(hillclimb(np.array([1,0.5]),f,dx,dy))