数学とか語学とか楽しいよね

フランス語、ドイツ語、ロシア語、アラビア語、オランダ語、英語、スペイン語、ラテン語とか数学とか数値計算(有限要素法、有限体積法、差分法、格子ボルツマン法、数理最適化、C++コード付き)とか勉強したことをまとめます。右のカテゴリーから興味のある記事を探してください。最近はクラシックの名演も紹介しています。Amazonアソシエイトを使用しています。

【アラビア語】アラビア語をはじめました!!!

はじめに

アラビア語をはじめました!!!実は第二外国語として学んだので厳密に言うと再開ですね。しかし、そのときはフラフラになりながら派生形をやって終わった、という程度です。単語も全然覚えていません。ただアラビア文字は全部読めるし書ける状態です。授業では『現代アラビア語入門』という参考書を使っていましたが、今見返してみても、初心者がこの本を勉強してアラビア語が読めるようになるとは到底思えません。一冊目としては解説が簡潔すぎます。なので、今回は別の参考書からはじめてみようと思います。ちなみに辞書はHans Wehrの"Arabic English Dictionary of Modern Written Arabic"でした。さすがにこれを初学者に与えるのはやりすぎでしたね。辞書も別のものを使っていきたいと思います。

現代アラビア語入門

現代アラビア語入門

Arabic English Dictionary of Modern Written Arabic

Arabic English Dictionary of Modern Written Arabic

目標

目標は、アラビア語Wikipediaを読めるようになることです。そのために、取りあえず文法は全体をやって、単語を確実に1000以上覚えたいです。欲を言えばコーランや本まで読めるようになりたいですね。

どんな本を読むか

文法

まず文法をやるということで、『ニューエクスプレス アラビア語』をひと通り読みました。100ページちょっとでアラビア語の文法や単語を派生形に入るかどうかまで解説してくれるのでとても便利です。下の記事も参考にしてください。

ニューエクスプレス アラビア語

ニューエクスプレス アラビア語


ここからは予定ですが、次に読む本としては『アラビア語の入門』ですね。この本は、『ニューエクスプレス アラビア語』と同じような範囲を扱っていますが、ページをたくさん使ってアラビア語に親しませようと努力してくれます。この続編が『ステップアップ アラビア語の入門』です。最近新版が出てうれしいです。派生形にかなり力を入れており、「この本を終えれば、アラビア語の初級は完成。」とのことです。また、副読本として『ニューエクスプレス アラビア語』と同じ著者による『アラビア語表現とことんトレーニング』がよさそうです。文法と単語を無理なく定着させるのに使えそうです。こちらも文法全体を扱っています。初級文法が固まったらまとめとして『アラビア語文法ハンドブック』を通読したいと思います。分厚いですが非常に見やすいレイアウトです。

アラビア語の入門 (<CD+テキスト>)

アラビア語の入門 ()

新版 ステップアップ アラビア語の入門

新版 ステップアップ アラビア語の入門

アラビア語表現とことんトレーニング

アラビア語表現とことんトレーニング

アラビア語文法ハンドブック

アラビア語文法ハンドブック

単語

文法を学びながらある程度単語を覚えてから単語帳を使おうと思います。最初から単語帳を使うと未知語ばかりでしんどいからです。まずは、知人が薦めてくれた『これなら覚えられる! アラビア語単語帳』を読もうと思います。基本単語約1400語が載っているようです。例文もあってよいですね。もう一冊は『例文で学ぶ アラビア語単語集』です。この本はこの前出たばかりで、本屋で見てみたところかなり有望です。こちらは2500語載っているようです。『これなら覚えられる! アラビア語単語帳』のあとに取り組むとちょうどよさそうです。例文があるのはやはりよいものです。最終的には"A Frequency Dictionary of Arabic: Core Vocabulary for Learners"まで行きたいですが、時間がかかるでしょう。頻度順にアラビア語単語5000語を並べた本です。例文もひとつずつ付いています。

NHK出版CDブック これなら覚えられる!  アラビア語単語帳

NHK出版CDブック これなら覚えられる! アラビア語単語帳

例文で学ぶ アラビア語単語集

例文で学ぶ アラビア語単語集

辞書

前回の反省から辞書は語根順ではない『パスポート 初級アラビア語辞典』を使います。文字でひけるし、用例も載っています。「見出し語は約4,200語。また最も頻度の高い基本単語として約500語を選び出し、窓見出しとなっている。」とのことなので、将来的に単語帳としても使えそうです。オンラインの辞書としては「アラジン」が便利です。アラビア語をクリックで入力できるのでとても楽です。用例もあり、派生語や派生形も出てくる、活用表も出てくる、などなどすさまじい機能があります。私は全然使いこなせていません。よくわからない単語はアラジンにぶちこむと大概わかります。下手な辞書より何倍も使い勝手がよいです。

パスポート 初級アラビア語辞典

パスポート 初級アラビア語辞典

おわりに

今後の見通しとしては、『ニューエクスプレス アラビア語』の復習(出てくる単語を全部覚えてしまいたい)と『アラビア語の入門』を読み進めることですね。『ステップアップ アラビア語の入門』までは失速する前に一気に行きたいところです。

【ドイツ語】『ドイツ語の小説を読む〈1〉ベル:きまぐれな客たち』【14冊目】

はじめに

佐藤清昭 著『ドイツ語の小説を読む〈1〉ベル:きまぐれな客たち』を読んだので、その紹介をしようと思います。ドイツ語で書かれた短編小説を精読することによって、読むためのドイツ語を身につけるための本です。

ドイツ語の小説を読む〈1〉ベル:きまぐれな客たち

ドイツ語の小説を読む〈1〉ベル:きまぐれな客たち

どんな本?

初級文法を学んだあとに、実際にドイツ語で書かれた短編小説を精読することによって、読むためのドイツ語を身につけるための本です。「ドイツ語の初級文法を一通り終え、「さて、次のステップは?」」と考えている人向けの本です。精読する短編小説は、Heinrich Böll(ハインリッヒ・ベル)の『きまぐれな客たち』("Unberechenbare Gäste")です。主人公の家に次から次へと珍客(動物たち)が訪れるというユーモラスな作品です。小説自体の長さは8ページ分で、それをじっくり精読していきます。ドイツ語学者である関口存男による、関口文法の影響を受けているのが本書の特徴です。この本自体のページ数は144ページでコンパクトです。

構成は?

小説全体を4部に分けて、さらに各部を約10課に分割しています。各課は見開き1ページで完結しており、読みやすいです。まず、小説の本文があります。10行満たないぐらいで、長すぎなくてよいです。その後丁寧な文法・構文・語彙の説明があります。ここを読めば独学でも十分に文章の構造を理解できます。そして、最後に日本語訳があります。私の大好きな古き良き構成です。また、ところどことに文法コラムが挿入されます。後でも述べるように、例文選びにセンスがあり、言語学的な説明も豊富で非常に読み甲斐があります。

f:id:mutsumunemitsutan:20190904204140j:plain:w500

良い点

  • 文法を使って実際のドイツ語をどうやって読むかが懇切丁寧に解説される

やはり理論だけわかっても読めないものです。このように実地で丁寧に鍛えてくれる本は貴重です。

  • 豊富な例文を用いて文法の解説をしてくれる

例文選びにセンスがあり、言語学的な説明も豊富です。ドイツ語に対して、外面(文法)と内面(文脈)からせまるアプローチが面白いです。このへんが関口存男の影響の一部かもしれません。

  • 話自体が面白い

ベルの文章はユーモラスで読んでいて楽しくなります。また、構文が複雑すぎないのもよいです。

  • あまり長くないため息切れせずに読み切ることができる

どんな本でもやはり読み切らないと効果は出ません。

  • 文脈の中で単語・熟語を覚えることができる

悪い点

  • 144ページと短いためすぐに読み終わってしまい人によっては物足りないかも

これを補うために下記のように姉妹編があります。

おわりに

非常に読んでいて爽やかで軽やかな本でした。ドイツ語で短編小説を読み終えたというのは自身に繋がります。次は姉妹編の、佐藤清昭 著『ドイツ語の小説を読む〈2〉ベル:ある若き王様の思い出』を読んでいこうと思います。

【ドイツ語】『文法復習やさしい独文解釈』【13冊目】

はじめに

有田潤 著『文法復習やさしい独文解釈』を読んだので、その紹介をしようと思います。初級と中級をつなぐ名著でした。

文法復習やさしい独文解釈

文法復習やさしい独文解釈

どんな本?

初級文法は学んだけれど、まだドイツ語の本は読めなくて次に何を勉強しようかと考えている人が読むべきドイツ語読解の参考書です。もっと詳しく言うと、「初級文法を一度学んだけれども、もう一度初級文法をやるのは面白くないし、かといって本を読むには単語や熟語の力が足りない人が中級へと到達すること」を想定したドイツ語読解の参考書です。取り上げられている題材は、芸術、地理、経済、歴史、情景描写など様々で、将来ドイツ語で本を読む時に役に立ちます。ページ数は188ページでコンパクトです。この本を読み切ると、中級ドイツ語の参考書が使えるようになります。例えば『独文解釈の研究』や『独文解釈の秘訣』などが読めるようになります。これらは「独学で大学教師よりドイツ語ができるようになる方法」でもおすすめされているドイツ語解釈の参考書です。是非読みましょう。私もまだ途中ですが…

独文解釈の研究

独文解釈の研究

独文解釈の秘訣―大学入試問題の徹底的研究 (1)

独文解釈の秘訣―大学入試問題の徹底的研究 (1)

独文解釈の秘訣―大学入試問題の徹底的研究 (2)

独文解釈の秘訣―大学入試問題の徹底的研究 (2)


構成は?

『独文解釈の研究』と非常に似ており、50課から構成されています。各課は、まず題材の文章(数百語程度の長さ)、単語・熟語の意味、構文・文法の解説、全訳、そして最後に作文となっています。各課ごとに注目する文法事項が示されていて重点的に解説してくれます。私はこの古き良き構成が大好きですし、読解力をつけるうえで一番効率的だと信じています。後ろにいくほど文章が難しくなっていきますが、難しくなりすぎることはありません。作文に関しては、私はもちろん答えのほうを見て読解として活用しています。

f:id:mutsumunemitsutan:20190823000037j:plain:w500

良い点

  • 単語・熟語を大量に覚えることができる

かなり丁寧な語釈がついているので辞書をひく手間が減らせます。

  • ドイツ語らしい構文に慣れることができる

まさにドイツ語らしい構文(関係代名詞や指示代名詞でどんどんつなぐ、接続詞を省略)のオンパレードで力がつきます。


  • 構文・文法の説明が丁寧

初級文法もしっかり復習できます。

  • 様々なタイプの文章が収録されている

ハイネ、モーツァルトベートーヴェンなどの芸術の話や、地理、経済、歴史、情景描写などが収録されており、将来ドイツ語で本を読む時に役に立ちます。もちろんそれぞれの話自体も興味深いです。

悪い点

悪い点は無いです!これは素晴らしい本ですね。

おわりに

実は著者の有田潤氏は『初級ラテン語入門』の著者でもあるのです!これはびっくりですね。ラテン語を研究するためにはドイツ語、フランス語は必須ということでしょうか。私も見習いたい限りです。

初級ラテン語入門

初級ラテン語入門

これでしっかり通して読み終わったので、一日一周の復習ができるようになりました。文章だけ読むと通しで1時間弱ですね。30回ぐらい読みたいです。

【最適化】準ニュートン法(BFGS公式とアルミホ条件による直線探索)C++コード付き

はじめに

今回は無制約最小化問題に対する数値解法(反復法)である、準ニュートン法(BFGS公式とアルミホ条件による直線探索)のC++コードを公開します。例題として2変数関数

 \min f(x_1,x_2) = \frac{1}{2} (x_1-x_2^2)^2 + \frac{1}{2} (x_2-2)^2

を考えます。最適解は  \boldsymbol{x}^*= (4, 2) です。反復法とは、適当な初期値  \boldsymbol{x}^{0} を定め、

 \boldsymbol{x}^{n+1} = \boldsymbol{x}^{n} + \alpha \boldsymbol{d}

という漸化式によって値を更新していき、最終的に最適解  \boldsymbol{x}^* へと収束させるような方法のことです。ここで、 \boldsymbol{d} が探索方向、 \alpha がステップのサイズです。つまり、反復法を構成しているのは探索方向とステップのサイズです。準ニュートン法の場合、ニュートン法とは異なりステップのサイズは直線探索を行います。

参考文献

参考文献は『最適化とその応用』pp.156-165、『最適化と変分法』pp.49-57、『最適化法』pp.119-122それと以下のPDFファイルです。今回は、特に『最適化とその応用』のp.158の例6が非常に助けになりました。準ニュートン法の1ステップを具体的な数値を出しながら説明してくれています。これを丁寧に追えば、自分のコードが正しいか判断できます。『最適化とその応用』は他の話題でも具体的な数値例が出てくるので大好きです。

http://www-optima.amp.i.kyoto-u.ac.jp/~nobuo/Ryukoku/2002/course7.pdf

工学基礎 最適化とその応用 (新・工科系の数学)

工学基礎 最適化とその応用 (新・工科系の数学)

基礎系 数学 最適化と変分法 (東京大学工学教程)

基礎系 数学 最適化と変分法 (東京大学工学教程)

最適化法 (工系数学講座 17)

最適化法 (工系数学講座 17)

準ニュートン法とは?

準ニュートン法とは、ニュートン法を改良した手法です。ニュートン法の説明は以下の記事を見て下さい。

ニュートン法はまとめると

 \nabla^2 f(\boldsymbol{x}) \boldsymbol{d} = -\nabla f(\boldsymbol{x})

という連立一次方程式を解いて  \boldsymbol{d} を求め

 \boldsymbol{x}^{n+1} = \boldsymbol{x}^{n} +  \boldsymbol{d}

のように値を更新して収束させる手法です。ただし、ニュートン法はヘシアン  \nabla^2 f(\boldsymbol{x}) が正定値であれば降下方向(目的関数の値が減少する方向)を与えるのですが、一般にはヘシアンは常には正定値になりません。そこで、常に正定値となるようなヘシアンを近似する行列  B をヘシアンの代わりに使おう、というのが準ニュートン法です。つまり、方向は

 B \boldsymbol{d} = -\nabla f(\boldsymbol{x})

を解いて決めます。今回は連立一次方程式の解法として直接法であるガウスの消去法を使います。以下にC++コードがあります。

方向の次はステップのサイズですが、これはアルミホの条件で決めます。以下の記事を見て下さい。

さて最後に残ったのがヘシアンの近似行列  B です。これは、BFGS公式と呼ばれる式を用いて更新していきます。Broyden、Fletcher、Goldfarb、Shannoらによって独立に提案されました。以下のような式です。

 B^{n+1} = B^{n} -\frac{B^{n}\boldsymbol{s}^{n} (B^{n}\boldsymbol{s}^{n})^{T}}{(\boldsymbol{s}^{n})^{T} B^{n} \boldsymbol{s}^{n}} + \frac{\boldsymbol{y}^{n} (\boldsymbol{y}^{n})^{T}}{(\boldsymbol{s}^{n})^{T} \boldsymbol{y}^{n}}

ここで、 \boldsymbol{s}^{n} = \boldsymbol{x}^{n+1} - \boldsymbol{x}^{n} \boldsymbol{y}^{n} = \nabla f(\boldsymbol{x}^{n+1}) - \nabla f(\boldsymbol{x}^{n}) です。大事なのは、この  B がセカント条件と正定値条件を満たすことです。これらは、ヘシアンが持っていて然るべき性質をあらわしています。詳細は『最適化法』pp.120-121を見て下さい。

準ニュートン法の詳細は上記の参考文献を見てください。

C++コード

#include <iostream>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <vector>

using namespace std;

//function//
inline double f(double x, double y)
{
	return 0.5*(x-y*y)*(x-y*y)+0.5*(y-2.0)*(y-2.0);
}

//gradient x//
inline double fx(double x, double y)
{
	return x-y*y;
}

//gradient y//
inline double fy(double x, double y)
{
	return -2.0*y*(x-y*y)+y-2.0;
}

//matix vector multiplication Ax=b
inline void MatVec(double A[],double x[],double b[], int size)
{
	for(int i=0;i<size;i++)
	{
		b[i]=0.0;
		for(int j=0;j<size;j++)
		{
			b[i]+=A[i*size+j]*x[j];
		}
	}
}

//inner product//
inline double ip(double a[], double b[], int size)
{
	double value=0.0;
	for(int i=0;i<size;i++) value+=a[i]*b[i];
	return value;
}

//Newton method, Gaussian elimination//
inline double Newton(int n, double b[], double B[], double x[])
{	
	double A[n*n];
	
	//copy//
	for(int i=0;i<n;i++)
	{
		for(int j=0;j<n;j++)
		{
			A[i*n+j]=B[i*n+j];
		}
	}
	
	//forward elimination//
	for(int i=0;i<n-1;i++)
	{
		for(int j=i+1;j<n;j++)
		{
			int m=A[j*n+i]/A[i*n+i];
			
			for(int k=i+1;k<n;k++)
			{
				A[j*n+k]=A[j*n+k]-m*A[i*n+k];
			}
			
			b[j]=b[j]-m*b[i];
			
			A[j*n+i]=0.;
		}
	}
	
	//backward substitution//
	x[n-1]=b[n-1]/A[(n-1)*n+n-1];
	
	for(int i=n-2;i>=0;i--)
	{
		for(int j=i+1;j<n;j++)
		{
			b[i]-=A[i*n+j]*x[j];
		}
		
		x[i]=b[i]/A[i*n+i];
	}
}

//dysdic, vector*vector=matrix, a*b=ab//
inline void dyadic(double a[], double b[], double ab[], int size)
{
	for(int i=0;i<size;i++)
	{
		for(int j=0;j<size;j++)
		{
			ab[i*size+j]=a[i]*b[j];
		}
	}
}

inline void Bupdate(double Bs[], double dgrad[], double B[], double sBs, double sdgrad, int size)
{
	double BsBs[size*size];
	double dgraddgrad[size*size];
	dyadic(Bs,Bs,BsBs,size);
	dyadic(dgrad,dgrad,dgraddgrad,size);
	
	for(int i=0;i<size;i++)
	{
		for(int j=0;j<size;j++)
		{
			B[i*size+j]=B[i*size+j]-BsBs[i*size+j]/sBs+dgraddgrad[i*size+j]/sdgrad;
		}
	}
}

//Armijo line search//
inline double Armijo(double x, double y, double d[])
{
	double alpha=1.0;//step size
	double rho=0.4;//contraction ratio
	double c1=0.001;//coefficient
	
	for(int i=0;i<10000;i++)
	{
		if(f(x+d[0]*alpha,y+d[1]*alpha)<=f(x,y)+c1*(fx(x,y)*d[0]+fy(x,y)*d[1])*alpha) break;
		else alpha=rho*alpha;
	}
	
	return alpha;
}

int main()
{
	int n=2;
	double err=pow(10.0,-10.0);//tolerance
	double x,y;//unknown values
	double xnew,ynew;//next time step values
	double d[n];//direction and step size
	double B[n*n];//approximate hessian
	double grad[n];//gradient
	double s[n];
	double dgrad[n];
	double alpha;//step size
	double Bs[n];
	double sdgrad;
	double sBs;
	
	//initial value//
	x=0.0;
	y=0.0;
	B[0*n+0]=1.0;
	B[0*n+1]=0.0;
	B[1*n+0]=0.0;
	B[1*n+1]=1.0;
	
	//iteration//
	for(int i=0;i<10000;i++)
	{
		cout<<i<<" "<<abs(fx(x,y))<<" "<<abs(fy(x,y))<<endl;
		
		//convergence check//
		if((abs(fx(x,y))<err)&&(abs(fy(x,y))<err)) break;
		
		grad[0]=-fx(x,y);
		grad[1]=-fy(x,y);
		
		//direction//
		Newton(n,grad,B,d);
		
		//step size//
		alpha=Armijo(x,y,d);
		
		//update for x//
		xnew=x+alpha*d[0];
		ynew=y+alpha*d[1];
		
		//update for B//
		s[0]=alpha*d[0];
		s[1]=alpha*d[1];
		dgrad[0]=fx(xnew,ynew)-fx(x,y);
		dgrad[1]=fy(xnew,ynew)-fy(x,y);
		MatVec(B,s,Bs,n);
		sdgrad=ip(s,dgrad,n);
		sBs=ip(s,Bs,n);
		
		Bupdate(Bs,dgrad,B,sBs,sdgrad,n);
		
		//system("pause");
		//update//
		x=xnew;
		y=ynew;
	}
	
	cout<<endl;
	cout<<x<<" "<<y<<endl;
	
    return 0;
}

【最適化】多変数ニュートン法 C++コード付き

はじめに

今回は無制約最小化問題に対する数値解法(反復法)である、多変数ニュートン法C++コードを公開します。例題として2変数関数

 \min f(x,y) = x^3 + y^3 -9xy +27

を考えます。最適解は  \boldsymbol{x}^*= (3, 3) です。反復法とは、適当な初期値  \boldsymbol{x}^{0} を定め、

 \boldsymbol{x}^{n+1} = \boldsymbol{x}^{n} + \alpha \boldsymbol{d}

という漸化式によって値を更新していき、最終的に最適解  \boldsymbol{x}^* へと収束させるような方法のことです。ここで、 \boldsymbol{d} が探索方向、 \alpha がステップのサイズです。つまり、反復法を構成しているのは探索方向とステップのサイズです。ニュートン法の場合、探索方向とステップのサイズは一度に求まります。お得ですね。

参考文献

参考文献は『これならわかる最適化数学』pp.88-93、『最適化と変分法』pp.44-49、『最適化とその応用』pp.151-155、それと以下のPDFファイルです。『これならわかる最適化数学』は特にビギナー向けで読みやすいです。
http://www-optima.amp.i.kyoto-u.ac.jp/~nobuo/Ryukoku/2002/course7.pdf

これなら分かる最適化数学―基礎原理から計算手法まで

これなら分かる最適化数学―基礎原理から計算手法まで

基礎系 数学 最適化と変分法 (東京大学工学教程)

基礎系 数学 最適化と変分法 (東京大学工学教程)

工学基礎 最適化とその応用 (新・工科系の数学)

工学基礎 最適化とその応用 (新・工科系の数学)

ニュートン法とは?

ニュートン法とは、目的関数(最小化しようとする対象)に対する1次の最適性条件

 \nabla f(\boldsymbol{x}) = \boldsymbol{0}

に対して、非線型方程式に対するニュートン法を適用したもの、と言えます。詳しく説明しましょう。現在  \boldsymbol{x}^{n} にいて  \boldsymbol{d} だけ進んだ時に1次の最適性条件を満たすとしましょう。つまり、

 \nabla f(\boldsymbol{x}+\boldsymbol{d}) = \boldsymbol{0}

です。我々が知りたいのは  \boldsymbol{d} です。上の式を1次までテイラー展開すると

 \nabla f(\boldsymbol{x}) + \nabla^2 f(\boldsymbol{x}) \boldsymbol{d}= \boldsymbol{0}

となります。 \nabla^2 f(\boldsymbol{x}) はヘシアンと呼ばれています。ここから  \boldsymbol{d} について解くと

 \boldsymbol{d} = -{\nabla^2 f(\boldsymbol{x})^{-1}} \nabla f(\boldsymbol{x})

となります。これが多変数ニュートン法における探索方向とステップサイズです。一度に求まります。多変数の場合は行列になるので逆行列をかける必要があります。これが1変数ニュートン法との違いです。この  \boldsymbol{d} を用いて

 \boldsymbol{x}^{n+1} = \boldsymbol{x}^{n} +  \boldsymbol{d}

のように値を更新して収束させればよいのです。ただし、上記のように逆行列を求めるのは計算コストが高いので、実際には

 \nabla^2 f(\boldsymbol{x}) \boldsymbol{d} = -\nabla f(\boldsymbol{x})

という連立一次方程式を解いて  \boldsymbol{d} を求めます。今回は連立一次方程式の解法として直接法であるガウスの消去法を使います。以下にC++コードがあります。

多変数ニュートン法の詳細は上記の参考文献を見てください。

C++コード

#include <iostream>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <vector>

using namespace std;

//function//
inline double f(double x, double y)
{
	return x*x*x+y*y*y-9.0*x*y+27.0;
}

//gradient x//
inline double fx(double x, double y)
{
	return 3.0*x*x-9.0*y;
}

//gradient y//
inline double fy(double x, double y)
{
	return 3.0*y*y-9.0*x;
}

//hessian xy//
inline double fxy(double x, double y)
{
	return -9.0;
}

//hessian fxx//
inline double fxx(double x, double y)
{
	return 6.0*x;
}

//hessian fyy//
inline double fyy(double x, double y)
{
	return 6.0*y;
}

//Newton method, Gaussian elimination//
inline double Newton(int n, double b[], double A[], double x[])
{	
	//forward elimination//
	for(int i=0;i<n-1;i++)
	{
		for(int j=i+1;j<n;j++)
		{
			int m=A[j*n+i]/A[i*n+i];
			
			for(int k=i+1;k<n;k++)
			{
				A[j*n+k]=A[j*n+k]-m*A[i*n+k];
			}
			
			b[j]=b[j]-m*b[i];
			
			A[j*n+i]=0.;
		}
	}
	
	//backward substitution//
	x[n-1]=b[n-1]/A[(n-1)*n+n-1];
	
	for(int i=n-2;i>=0;i--)
	{
		for(int j=i+1;j<n;j++)
		{
			b[i]-=A[i*n+j]*x[j];
		}
		
		x[i]=b[i]/A[i*n+i];
	}
}

int main()
{
	int n=2;
	double err=pow(10.0,-10.0);//tolerance
	double x,y;//unknown
	double d[n];//direction and step size
	double Hess[n*n];//hessian
	double grad[n];//gradient;
	
	//initial value//
	x=5.0;
	y=7.0;
	
	//iteration//
	for(int i=0;i<10000;i++)
	{
		cout<<i<<" "<<abs(fx(x,y))<<" "<<abs(fy(x,y))<<endl;
		
		//convergence check//
		if((abs(fx(x,y))<err)&&(abs(fy(x,y))<err)) break;
		
		grad[0]=-fx(x,y);
		grad[1]=-fy(x,y);
		
		Hess[0*n+0]=fxx(x,y);
		Hess[0*n+1]=fxy(x,y);
		Hess[1*n+0]=fxy(x,y);
		Hess[1*n+1]=fyy(x,y);
		
		Newton(n,grad,Hess,d);
		
		//update//
		x=x+d[0];
		y=y+d[1];
	}
	
	cout<<endl;
	cout<<x<<" "<<y<<endl;
	
    return 0;
}

【最適化】1変数ニュートン法 C++コード付き

はじめに

今回は無制約最小化問題に対する数値解法(反復法)である、1変数ニュートン法C++コードを公開します。例題として

 \min (x-2)^2

を考えます。もちろん、最適解は  x^*=2 です。反復法とは、適当な初期値  x^{0} を定め、

 x^{n+1} = x^{n} + \alpha d

という漸化式によって値を更新していき、最終的に最適解  x^* へと収束させるような方法のことです。ここで、 d が探索方向、 \alpha がステップのサイズです。つまり、反復法を構成しているのは探索方向とステップのサイズです。ニュートン法の場合、探索方向とステップのサイズは一度に求まります。お得ですね。

参考文献

参考文献は『最適化と変分法』pp.44-49、『最適化とその応用』pp.151-155、それと以下のPDFファイルです。どちらの本も非常に参考になります。いつも読んでいます。
http://www-optima.amp.i.kyoto-u.ac.jp/~nobuo/Ryukoku/2002/course7.pdf

基礎系 数学 最適化と変分法 (東京大学工学教程)

基礎系 数学 最適化と変分法 (東京大学工学教程)

工学基礎 最適化とその応用 (新・工科系の数学)

工学基礎 最適化とその応用 (新・工科系の数学)

ニュートン法とは?

ニュートン法とは、目的関数(最小化しようとする対象)に対する1次の最適性条件

 \nabla f(x) = 0

に対して、非線型方程式に対するニュートン法を適用したもの、と言えます。詳しく説明しましょう。現在  x^{n} にいて  d だけ進んだ時に1次の最適性条件を満たすとしましょう。つまり、

 \nabla f(x+d) = 0

です。我々が知りたいのは  d です。上の式を1次までテイラー展開すると

 \nabla f(x) + \nabla^2 f(x) d= 0

となります。 \nabla^2 f(x) はヘシアンと呼ばれています。ここから  d について解くと

 d = -\frac{\nabla f(x)}{\nabla^2 f(x)}

となります。これが1変数ニュートン法における探索方向とステップサイズです。一度に求まります。多変数の場合は行列になるので逆行列をかける必要があります。この  d を用いて

 x^{n+1} = x^{n} +  d

のように値を更新して収束させればよいのです。詳細は上記の参考文献を見てください。

C++コード

#include <iostream>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <vector>

using namespace std;

//function//
inline double f(double x)
{
	return (x-2.0)*(x-2.0);
}

//gradient//
inline double fx(double x)
{
	return 2.0*(x-2.0);
}

//hessian//
inline double fxx(double x)
{
	return 2.0;
}

//Newton method//
inline double Newton(double x)
{	
	return -fx(x)/fxx(x);
}

int main()
{
	double err=pow(10.0,-10.0);//tolerance
	double x=10.0;//initial value
	double d;//direction and step size
	
	
	//iteration//
	for(int i=0;i<10000;i++)
	{
		cout<<i<<" "<<abs(fx(x))<<endl;
		
		//convergence check//
		if(abs(fx(x))<err) break;
		
		d=Newton(x);
		
		//update//
		x=x+d;
	}
	
	cout<<endl;
	cout<<x<<endl;
	
    return 0;
}

【最適化】最急降下法(アルミホ条件による直線探索)C++コード付き

はじめに

今回は無制約最小化問題に対する数値解法(反復法)である、最急降下法(アルミホ条件による直線探索)のC++コードを公開します。例題として

 \min (x-2)^2

を考えます。もちろん、最適解は  x^*=2 です。反復法とは、適当な初期値  x^{0} を定め、

 x^{n+1} = x^{n} + \alpha d

という漸化式によって値を更新していき、最終的に最適解  x^* へと収束させるような方法のことです。ここで、 d が探索方向、 \alpha がステップのサイズです。つまり、反復法を構成しているのは探索方向とステップのサイズです。今回は、探索方向としては最急降下法を、ステップのサイズとしてはアルミホ(Armijo)条件による直線探索を行います。

参考文献

参考文献は『最適化と変分法』pp.35-43、『最適化とその応用』pp.130-139、それと以下のPDFファイルです。どちらの本も非常に参考になります。いつも読んでいます。
http://www-optima.amp.i.kyoto-u.ac.jp/~nobuo/Ryukoku/2002/course7.pdf

基礎系 数学 最適化と変分法 (東京大学工学教程)

基礎系 数学 最適化と変分法 (東京大学工学教程)

工学基礎 最適化とその応用 (新・工科系の数学)

工学基礎 最適化とその応用 (新・工科系の数学)

最急降下法とは?

最急降下法とは、勾配が最もきつい方向、すなわちgradのマイナス方向に進んでいく手法のことを言います。勾配のきつい方向に進んでいくと底の部分にたどり着けるというわけです。ただし、場合によってはこの方法は効率がよくありません。そのような場合にも適用できる手法として共役勾配法などがあります。詳細は上記の参考文献を見てください。

アルミホ条件による直線探索とは?

アルミホ(Armijo)条件による直線探索とは、スタート地点での傾きに1より小さい係数をかけた傾きを使って直線をひき、その直線よりも関数の値が下回るようなステップサイズを採用する、というものです。上記のPDFの図がわかりやすいです。詳細は上記の参考文献を見てください。

C++コード

#include <iostream>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <vector>

using namespace std;

//function//
inline double f(double x)
{
	return (x-2.0)*(x-2.0);
}

//gradient//
inline double fx(double x)
{
	return 2.0*(x-2.0);
}

//steepest gradient//
inline double SteepGrad(double x)
{
	return -fx(x);
}

//Armijo line search//
inline double Armijo(double x, double d, double alpha, double c1, double rho)
{
	alpha=1.0;
	for(int i=0;i<10000;i++)
	{
		if(f(x+d*alpha)<=f(x)+c1*fx(x)*d*alpha) break;
		else alpha=rho*alpha;
	}
	
	return alpha;
}

int main()
{
	double err=pow(10.0,-10.0);//tolerance
	double x=10.0;//initial value
	double alpha;//step size
	double rho=0.4;//contraction ratio
	double c1=0.001;//coefficient
	double d;//direction
	
	
	//iteration//
	for(int i=0;i<10000;i++)
	{
		cout<<i<<" "<<abs(fx(x))<<endl;
		
		//convergence check//
		if(abs(fx(x))<err) break;
		
		//step size search//
		d=SteepGrad(x);//direction
		alpha=Armijo(x,d,alpha,c1,rho);//step size;
		
		//update//
		x=x+alpha*d;
	}
	
	cout<<endl;
	cout<<x<<endl;
	
    return 0;
}