【確率微分方程式】確率微分方程式を数値計算してみよう C++コード付き

はじめに

今回は確率微分方程式を数値的に計算するC++コードを紹介しようと思います。パスがランダムに進展していく様は見ていて非常に興味深いです。

確率微分方程式とは？

一般論

のような式のことを伊藤の確率微分方程式（Stochastic Differential Equation、SDE）と呼んでいます。ここで、 が時刻 における未知数（例えば株価や生物の個体数など）、 は時間、 はドリフト係数、 は拡散係数（ボラティリティ）、 は1次元の標準ブラウン運動で、この項が確率論的なゆらぎをあらわしています。標準ブラウン運動は、平均0、分散tの正規分布に従います。この事実を使うと、確率微分方程式を数値的に解くことができるようになります。

幾何ブラウン運動

今回扱う具体的な確率微分方程式は

としましょう。ここで、、 と設定しています。このようにドリフトと拡散係数が線型で を含むものを幾何ブラウン運動と呼びます。ウィナー過程とも呼ばれています。幾何ブラウン運動は数理ファイナンスや金融工学の分野で頻出します。例えば、ブラック・ショールズ方程式で使われています。幾何ブラウン運動の利点は線型で比較的解析が簡単であることと解析解が求まることです。上記の幾何ブラウン運動の解析解は

とあらわすことができます。ここで、 は初期値です。

オイラー・丸山スキーム

オイラー・丸山スキームは、確率微分方程式に対するもっともシンプルな数値計算手法です。常微分方程式におけるオイラー法と対応しています。収束次数は1/2次となかなかしょっぱいです。オイラー・丸山スキームは確率微分方程式を

のように近似します。ここで、 はタイムステップ、 は時間の刻み幅、 です。 は標準正規分布に従い、平均0、分散1を満たすものとしています。つまりブラウン運動の増分は、平均0、分散 の正規分布に従います。ここがとても大切です。数値計算では 、すなわち標準正規分布を発生させればよいのです。

計算結果

パラメータとしては、、、、、計算終了時刻は1としました。乱数はメルセンヌツイスターを使っています。計算結果を見ていきましょう。数値解と解析解の両方をプロットしています。青が数値解で、赤が解析解です。解析解は数値解と同じブラウン運動を使って計算します。

このように、ブラウン運動において、パスはギザギザしており微分不可能ですが、確率1で連続となります。値が増えたり減ったりランダムな動きをしていることがわかると思います。赤と青の線はほぼ重なっており、ちゃんと確率微分方程式が解けていることがわかりました。シード値を変えることにより異なるパスが実現します。いろいろ変えて遊んでみてください。

C++コード

#include <random>
#include <iostream>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <string>
#include <sstream>

using namespace std;

inline double drift(double x, double t, double mu)
{
    return x*mu;
}

inline double diffusion(double x, double t, double sigma)
{
    return x*sigma;
}

int main()
{
    int seed=80;
	mt19937 mt(seed);
    normal_distribution<> gauss(0.0, 1.0);

    ofstream file;
    file.open("path.txt");

    double mu=1.0;
    double sigma=sqrt(2.0);
    double x=1.0;
    double dt=0.0001;
    double TotalTime=1.0;
    double t;
    double exact=1.0;
    double random;
    double W=0.0;

    file<<0.0<<" "<<x<<" "<<exact<<endl;
    
	for(int i=0;dt*double(i)<TotalTime;++i)
	{
        t=dt*double(i);
        random=gauss(mt);
        W+=random*sqrt(dt);
        x=x+drift(x,t,mu)*dt+diffusion(x,t,sigma)*sqrt(dt)*random;
        exact=exp(sigma*W);

	    file<<t<<" "<<x<<" "<<exact<<endl;
	}

    return 0;
}

参考文献

『微分方程式による計算科学入門』の第4章を参考にしています。この本は確率微分方程式の数値解法に関して解説してある数少ない和書で重宝しています。

微分方程式による計算科学入門

あとは毛利光希氏によるこちらの資料です。
http://www-mmc.es.hokudai.ac.jp/~masakazu/Software/Manual/SDE_Solver/Manual.pdf

幾何ブラウン運動の解析解の導出は『確率微分方程式入門』のpp.77-78の例題4.3に載っています。伊藤の公式を適用するだけです。

確率微分方程式入門 ―数理ファイナンスへの応用― (数学のかんどころ 26)

【アラビア語】『アラビア語の入門』【2冊目】

はじめに

『ニューエクスプレス　アラビア語』の次の教材として、本田孝一 著『アラビア語の入門』を読んだのでその紹介をしようと思います。派生形の手前までの文法を丁寧に解説してくれる本です。

アラビア語の入門 ()

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

どんな本？

まったくの初心者がアラビア語をはじめるのに最適な本です。文法と会話のバランスがよく、実地で使うことを想定しています。アラビア文字の解説も丁寧で読み方、書き方を習得できます。文字自体の書き方やつなげかたまで丁寧です。扱っている範囲は派生形の手前ぐらいまでです。初心者がいきなり派生形を勉強しても挫折してしまうのでこれでちょうどよいと思います。そこまで初心者を挫折することなく導こうとする気概が感じられます。だいたい『ニューエクスプレス　アラビア語』と同じような範囲です。ページ数は202ページでコンパクトです。著者の本田氏はアラビア書道で有名な方です。なので、表紙やいたるところに氏の作品を見ることができます。

構成は？

全部で20課から構成されています。各課のページ数はまちまちですが10ページぐらいです。最初にキーフレーズや大事な文法事項が説明されて、それをテキスト（会話）や例文を通して学んでいくスタイルです。はじめて出てくる単語の意味もかいてあります。

良い点

• コンパクトで十分に読み切ることができる

読み切ることは大事です。さらに、図が多く、レイアウトも余裕があり目が疲れません。アラビア文字も十分大きく読みやすいです。

• 文法を小出しにしてくれる

文法を少しづつ小出しにしてくれるので、そのおかげで頭がパンクすることなく読み進められます。例えば、動詞はかなり後半まで読まないとでてきません。表に全ての活用を与えて覚えろ、ということはないです。もちろん文法自体の説明もスマートでわかりやすいです。

• 覚えるべき単語のリストがついている

基本単語（名詞と形容詞、名詞は複数形付き）が200、基本動詞が100のリストが載っています。初心者のうちは何を覚えたらよいかわからないのでとても重宝しています。全部ちゃんと覚えようと思います。

悪い点

• アラビア語初級文法すべてが載っているわけではない

入門としてはこのほうが負担が少なくて私はよいと思います。次は『ステップアップ　アラビア語の入門』でしっかり派生形を勉強してアラビア語初級文法の完成を目指しましょう！

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

• ある程度いろいろな語彙が出てくる

いろいろな語彙に触れられるのはよいことですが、少々多めかもしれません。なので、一周目で覚えきろうとせずに何度も読むとよいと思います。

おわりに

『ニューエクスプレス　アラビア語』を読んだ後にやったので非常にスムーズに読めました。まずは『ニューエクスプレス　アラビア語』をやるのがよさそうです。次はいよいよ『ステップアップ　アラビア語の入門』です！余裕があれば『アラビア語表現とことんトレーニング』を並行して読み進めたいです。

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

【アラビア語】『ニューエクスプレス アラビア語』【1冊目】

はじめに

最近アラビア語をはじめました。まずは文法を勉強しようと思い、竹田敏之 著『ニューエクスプレス　アラビア語』を読んだのでその紹介をしようと思います。まったくの初心者がアラビア語をはじめるのに最適な本です。

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

どんな本？

まったくの初心者がアラビア語をはじめるのに最適な本です。文法と会話のバランスがよく、実地で使うことを想定しています。アラビア文字の解説も丁寧で読み方、書き方を習得できます。扱っている範囲は派生形の手前ぐらいまでです。しかし、初心者がいきなり派生形を勉強しても挫折してしまうのでこれでちょうどよいと思います。だいたい『アラビア語の入門』と同じような範囲です。ページ数は153ページでコンパクトです。

アラビア語の入門 ()

構成は？

全部で20課から構成されています。各課は4ページで、最初の見開き1ページ目でテキスト（会話）が提示されます。左ページにアラビア語、右ページに和訳が載っており、対訳のようになっています。下の方には新しく出てきた単語の意味と慣用表現が載っています。これはなかなか便利です。一々辞書をひいているとどうしても時間がかかって勉強意欲が減退するからです。次のページの見開きは文法説明にあてられています。限られたスペースの中で例文を示しつつ非常にわかりやすく解説してくれます。2課ごとに練習問題もついています。しかし、私は練習問題はやらないようにしています（特に1周目）。面倒くさくなってやる気がなくなるからです。勉強を続けるために徹底的に楽をしましょう。また、単語力アップ、表現力アップというページも時折あり、語彙や表現を強化できます。巻末にはこの本に出てくるアラビア語の単語が集めてあります。便利なのは活用した形（私は～する、のような形）がそのまま載っていることです。これは初心者にはうれしいです。

良い点

• コンパクトで十分に読み切ることができる

読み切ることは大事です。

• 文法の説明がスマートでわかりやすい

初心者には十分です。さらに勉強するには『アラビア語文法ハンドブック』あたりでしょうか。

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

• ある程度語彙をしぼってあるのでパンクしないですむ

容赦のない語学書を読んでいると、最初から小難しい語彙がわんさか出てきて本質でない部分で挫折することになってしまいます。初級のうちは語彙をしぼり、文法の概観と骨格を学ぶことが大事だと思います。

悪い点

• アラビア語初級文法すべてが載っているわけではない

入門としてはこのほうが負担が少なくて私はよいと思います。この本にうまく接続できそうなのが『ステップアップ　アラビア語の入門』です。

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

• テキスト（会話）が無難で面白くない

はっきりいって普通です。パスポートを途中で無くすぐらいです。せっかく20課あるので何かストーリーが欲しかったですね。例えば『ドイツ語、もっと先へ！』はストーリーがぶっとんでいました。落語したり言語学を語りだしたり。

ドイツ語、もっと先へ！

おわりに

スムーズにアラビア語に入門できました。最初からこの本で勉強したかったですね。次は、『アラビア語の入門』を読みたいと思います。『ステップアップ　アラビア語の入門』まで頑張ります！

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

はじめに

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

現代アラビア語入門

Arabic English Dictionary of Modern Written Arabic

目標

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

どんな本を読むか

文法

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

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

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

アラビア語の入門 ()

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

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

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

単語

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

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

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

辞書

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

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

おわりに

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

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

はじめに

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

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

どんな本？

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

構成は？

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

良い点

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

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

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

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

• 話自体が面白い

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

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

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

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

悪い点

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

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

おわりに

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

ドイツ語の小説を読む〈2〉ベル:ある若き王様の思い出

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

はじめに

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

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

どんな本？

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

独文解釈の研究

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

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

構成は？

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

良い点

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

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

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

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

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

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

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

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

悪い点

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

おわりに

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

初級ラテン語入門

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

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

はじめに

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

を考えます。最適解は です。反復法とは、適当な初期値 を定め、

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

参考文献

参考文献は『最適化とその応用』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)

準ニュートン法とは？

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

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

という連立一次方程式を解いて を求め

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

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

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

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

ここで、、 です。大事なのは、この がセカント条件と正定値条件を満たすことです。これらは、ヘシアンが持っていて然るべき性質をあらわしています。詳細は『最適化法』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;
}