島田商業高等学校講演レジメ

  1. 静岡理工科大学情報学部の紹介
  2. 本日のプレゼン資料(PDFファイル)
  3. HTML5の機能
  4. フォームとPHPの連携
  5. Webデザイン特別プログラム製作物
  6. Webアプリケーションの例
    1. シューティングゲーム(リンク無)
    2. 名刺発行システム(リンク無)
    3. 関数グラフ&数表出力
    4. 電気料金節約シミュレーター

[模擬講義] コンピュータの計算は正しくない!

オープンキャンパス模擬講義

「コンピュータの計算は正しくない!」

静岡理工科大学 情報学部
コンピュータシステム学科
幸谷 智紀(こうや とものり)


1.「数」とは?

 現在のICT社会で使用される「データ(data)」は例外なく全て「数」を表現したものです。小学校では自然数,分数(有理数),小数を,中学校では負の数と整数を,高校では無理数を習いますが,それらの数を総称して「実数(real number)」と呼びます。実数を幾何学的に表現する時には直線=数直線を使用します。全ての数は,この数直線上のどこかに「点(point)」として存在しています。

 自然数,整数,有理数は全て分数として表現できる数です。下記の図に示す通り,分数は無数に存在します。

 ちょっと脱線しますが,下記のような問題はご存知でしょうか?

 解答については検索して調べて頂ければ何らかの手かがりが得られるでしょうが,実際の分数の計算はこのようにしてはいけませんね。どのように計算できるのか,ちょっと試してみましょう。

分数の計算フォーム

 実数にはもう一つ,分数としては表現できない「無理数(irrational number)」があります。有名なものとしては円周率\pi = 3.1415...や平方根\sqrt{2} = 1.4142...があります。これも実数の一種なので,数直線上の点として実在しています。

 有理数(分数)を小数で表現すると,次のように必ず同じ数のパターンが繰り返し出てきます。

\[\begin{split}
\frac{1}{2} &= 0.50000\cdots = 0.5\dot{0} \\
\frac{1}{3} &= 0.33333\cdots = 0.\dot{3}
\end{split}\]

 無理数では同じパターンの数の羅列が永久に続くことはありません。ちょっと計算してみましょう。

Try! MPFR

 従って,実数を一言で表現すると「循環したりしなかったりする)無限桁の小数として表わすことのできる数」ということになります。何故無限桁が必要になるのかは,微分積分(数学II, III)に登場する「極限」というものと深い関係がありますが,それは大学に入ってからじっくり考えてみて下さい。

 現在の科学技術では実数を土台とする計算が不可欠です。そしてほとんどの計算はコンピューターの中で行います。では,次にコンピューターの中身について見ていくことにしましょう。

2.コンピュータの仕組みと「データ」

 今のコンピューターの筐体とぱかっと開けると,マザーボードと呼ばれるデジタル回路基板の上に様々な集積回路が載っていることが分かります。下記は本研究室で故障したマシンのマザーボードの写真です。

 マザーボードの上に乗っている集積回路のうち,重要なものがCPU(Central Processing Unit)とメモリ(RAM)です。コンピューターで行われる処理は,全てのデータが「数」として表現されてメモリに格納され,処理内容に応じてメモリ内のデータが読み出されて処理され,またメモリに書き戻されます。

 このCPUとメモリとのやり取りに要する時間,CPU内の処理(計算)に要する時間が短ければ短いほど,コンピューターの性能が高いと言えます。とても高速な処理が可能になっていますが,困ったことに,メモリの大きさは有限で,処理時間も短いとはいえゼロにはなりません。「無限の長さを持つ実数」をそのままコンピューターのメモリに格納することは不可能ですし,仮にできたとしても,無限の長さの処理時間を要することになります。「無限大(∞)」とは果てのないことを意味し,それ故に完ぺきな数学,特に微分積分のような解析学という体系が成り立つわけですが,これではコンピューターには実数の計算が出来ないことになってしまいます。

3.「誤差」を含まざるを得ない科学技術計算

 仕方がないので,コンピューターの中では無限桁の小数を適当な長さに切って,有限桁の小数(=有理数の一種)として扱うことにします。この適当長さに切る操作のことを「丸め(Round-off)」と呼び,丸めに伴って起きる元の正しい値(真値)とのズレを「誤差(error)」と呼びます。

 この丸めに伴う誤差を,関数グラフアプリを使ってみていくことにしましょう。

MPFRgraph

コンピューターは計算が高速なので,いくらでも長い桁の計算ができそうですが,長くなればなるほど計算時間はかかるようになりますし,科学技術開発競争を世界中で行っている昨今,あまり誤差ばかり気にして長い桁の計算ばかりすることは難しいのが現状です。従って,最近のAIの中核技術であるDeep Learningではなるべく短い桁の小数で計算しようとしています。

 しかし,いくら計算が速くても,途中で誤差が拡大されて,不正確な値しか出てこないようでは困ります。一般には,丸める桁数を長くすれば誤差も小さくできるわけですが,とてもたちの悪い問題があることは良く知られていて,いくら桁を長く取っても足りなくなるということもあります。

5.誤差評価付き計算の例

 性質の良くない問題の例として,「複雑系(Chaos system)」というものがあります。一番シンプルなものとして,ロジスティック写像で定義される次のような数列x_0, x_1, ..., x_i, x_{i+1}, ...の生成問題を取り上げましょう。

 出発値x_0 = 0.7501として,次の式の右辺を計算して次のx_1を導出します。

\[ x_{i+1} := 4x_i (1 – x_i) \]

これで実数列\{x_i\}^{100}_{i=1}を,下記のCプログラム

を用いてIEEE754単精度(約7桁)及び倍精度(約15桁)で計算した結果は次のようになります。

同じ値になる名図ですが,x_{20}あたりでもう単精度の値は最初の一桁目しか合っていません。x_{30}以降はもう全く違う数です。

 このような場合はもっと桁数を増やして計算する必要があります。50桁計算すると下記のようになります。

100桁計算すると,

となり,大分同じ数に近づいてきたことが分かります。

 このように,丸めによる誤差の影響が大きく出る「悪条件問題」というものが,科学技術計算にはたまに出てきます。計算途中で誤差がどの程度は行っているのかを自動的に調べる「区間演算(interval arithmetic)」という技術もありますが,このような悪条件問題の場合は大した役には立たず,結局桁数を多くしないとまともな値を得ることはできません。下記の例は,MPFIと呼ばれる多倍長計算を区間演算に利用できるソフトウェアを使ったロジスティック写像の掲載例です。一応,プログラムもつけておきましょう。

6.まとめ

 以上,まとめますと

  1. 実数は無限桁の小数として表現される数である。
  2. コンピューターのメモリは有限なので,無限桁の実数を格納することはできず,丸めて短い有限桁に収めなければならない。
  3. 丸めに伴う誤差が甚大な影響を及ぼす「悪条件問題」というものが存在する。
  4. 悪条件問題を解決するには桁を増やして対処するしかない。しかしこれは計算時間を要する解決策である。

となります。「高性能計算研究室」ではこのように長い桁の計算,多倍長計算を用いて様々な悪条件問題の解決に取り組んでいます。

星陵高等学校講演:オープンデータとWebアプリケーション

「オープンデータとWebアプリケーション」
静岡理工科大学 情報学部 コンピュータシステム学科
幸谷 智紀(こうや とものり)

0.【宣伝】静岡理工科大学情報学部について
静岡理工科大学情報学部のWebページ

1.オープンデータとは?
 「オープンデータ(Open Data)」とは,誰でも閲覧できる状態で広く公開されているデータのことです。今の時代はデータを扱う主体はコンピューターですので,オープンデータと言っても紙媒体ではなく,電子化されているもの,即ち,コンピューターが直接読み書きできるファイルの形式になっているものを指すことが殆どです。
 電子化されたデータは,コンピューター同士が繋がっているインターネット(The Internet)の上でやり取りしやすい形式ですので,公開するにしても,インターネット上で誰しも読める状態になっている方が,読む方としても都合が良いわけです。インターネットでは様々なサービスが展開されていますが,現在ではWeb(World Wide Web)を通じたサービスが普通になっていますので,オープンデータもWeb上で,つまり,データを持っている会社,行政組織,個人,各種団体のWebページ(ホームページ)上に置き,誰でもブラウザを通じて読むことができるようにしています。こうしたオープンデータを活用したWeb上のサービス,いわゆるWebアプリケーションが様々な形で展開されています。

 以下,オープンデータとは,Web上で誰でもブラウザから閲覧できるようになっている電子ファイル状態のデータである,と定義します。最近では,行政機関が積極的にデータを公開するようになっています。静岡県も積極的にオープンデータの活用を後押ししていますし,静岡県内の市,例えば富士宮市,静岡市,浜松市では下記のところでオープンデータにアクセスできるようになっています。

 オープンデータは,単なる公開情報という存在を超えるものです。コンピューターで利用しやすいファイルの形式になっているということは,そのままソフトウェアで加工できることになります。これをうまく活用することで,より便利なWebアプリケーションを作ることができます。

2.オープンデータをどのように活用するか?
 ここでは,簡単なオープンデータの活用方法を具体例で示します。

 皆さんは,「国民の祝日」,例えば正月,成人の日,建国記念日,春分の日,がどこでどのように決まっているか,ご存知ですか? 検索すればすぐに分かりますが,これは政府が法律に基づいて決めているもので,内閣府に解説ページがあります。ただ,これは情報を人間に見やすく表示しているだけで,余り使い勝手の良いものではありません。そこで,近年の国民の祝日を表計算ソフトで扱いやすい形式で公開しています。これが国民の祝日のオープンデータです。この形式であれば,表計算ソフトだけでなく,様々なソフトウェアに利用することができるようになります。例えば,下記の図は,FullcalendarというJavaScriptで構築されたスケジュール帳に,国民の祝日を上のURLからファイルを読み取り,自動的に挿入したものです(サンプル)。

 但し,これはあまり使い勝手の良いものではありません。元々の国民の祝日ファイルは,単なるテキストファイルなので,データそのものを使いやすい形で形成し,流し込む処理が必要になります。例えば,上記のFullcalendarに取り込むためには,次のようなPHPスクリプト(プログラム)の手助けが必要です。

 このように,直接データを取得して自動的に加工できるようになると,人間を介しての作業が不要になります。

 オープンデータの利用のためには,もう少し,利用者の立場に立って,使い勝手の良い形式で必要な分量のデータを受け取れるようにする必要があります。最悪な形式は,いわゆる「ネ申Excel(三重大学・奥村晴彦)」と呼ばれる,書類をそのまま表計算ソフトウェアで成形しただけのシロモノで,人間が手作業でチェックするには便利かもしれませんが,コンピューターにとっては最も加工の難しいものです。オープンデータは人間よりも,コンピューターにとって使い勝手の良いシンプルな形式であることが望まれます。

 しかし,もっと使い勝手の良いオープンデータは,情報を加工したい側が欲しいデータだけをいつでも取得できるようになっているものです。そのためには,上記のFullcalendarに祝日データを渡すプログラムのように,オープンデータと利用者の間で仲介し,利用者の要求に応じた返答を返すような仕組みを作らなくてはなりません。これを実現しているのが,Web API (Application Programming Interface)と呼ばれるものです。

3.Web APIの実例
 最も良く知られたweb APIは,GoogleやYahoo! Japan,Bingといった検索エンジンです。

 どの検索エンジンも,URLの中に検索すべき単語を埋め込んであります。例えばGoogleの場合は

https://www.google.co.jp/#q=星陵高等学校

となっています。このように,ブラウザからアクセスする際に,利用者が欲しいデータの範囲(この場合は「星陵高等学校」を含むWebページの一覧)を指定することで,必要のないデータまで取り込むことがなくなります。

 現在では様々なデータを,Web APIを通じて取り込むことができるようになっています。Google map(地図情報)やぐるなび(グルメ情報)なども,細かくデータを指定して取り込むことができるようになっています。

 その他,下記のようなWeb APIもあります。

4.openBDを基盤とする書誌データの活用
 オープンデータとWeb API活用の事例として,現在本研究室で取り組んでいるopenBDについて紹介します。
 2017年初めに広く公開された本のオープンデータ化を進めるプロジェクトで,2017年6月現在では約90万件の本のデータが取得できます。今のところ,本のISBN番号一覧と,ISBN番号に紐づけされた書誌データ(表紙画像も含む)を取得できるWeb APIが提供されており,基本的にはユーザーがダウンロードして使用する全データを提供することを主目的としています。
 本研究室でも,簡単なWebのインターフェースを提供しています。

openBD API

本一冊一冊ごとに異なるISBN番号が振られており,それが分かればその本のデータだけを入手することが可能です。しかし,普通本を買うときにISBN番号で指定することはありません。普通は,本のタイトル,著者名,出版社名,定価等で探すものです。

 それを可能にするため,本研究室では全件検索ができるような仕組みを作り,公開しています。使い勝手はまだ改良の余地がありますが,これも改良されたopenBDのWeb API化の一つと言える仕事になっています。

5.おわりに
 インターネットが当たり前になった今では,様々なデータを自由に取得でき,もっと便利に使うためのWeb APIの利用が盛んになっています。今後は,これらのオープンデータやWeb APIを通じてより高度なデータ解析を自動的に,例えばAIを使って行うことが普通になっていきます。
 人間の側はもっと賢く,これらのオープンデータやWeb API,そしてその分析結果を使って高度な知を使いこなしていくことになります。プログラミングを行うだけでなく,それをより良い方向に生かす技術を,今後の大学生活で身につけて下さい。

大学見学会「高性能計算研究室」

高性能計算とWebアプリケーション

静岡理工科大学 コンピュータシステム学科 幸谷こうや 智紀とものり

1.概要

 データの入出力をWeb経由で行うためのソフトウェアの開発,すなわち「Webプログラミング」の事例をいくつか紹介します。関数の計算式を入力することでグラフや数表を自動的に出力するアプリケーションや特別プログラムの製作物を体験して頂く予定です。

2.「高性能計算」研究室とは?

 現在の科学技術計算に欠かすことのできないのが,コンピュータ上で現象を再現するコンピュータシミュレーション技術です。このコンピュータシミュレーションの信頼性を高め,かつ,効率的に実行するためには,計算のための理論的枠組み,計算を効率よく行うためのアルゴリズム,そして堅牢かつ高性能なコンピュータが欠かせません。このようなコンピュータシミュレーションを下支えする学問分野を「高性能計算(High  Performance Computing, HPC)」と呼びます。詳細はこちらの記事をご覧ください。

 本学の「高性能計算研究室」はその中でも特に,ソフトウェアを使ってシミュレーションの信頼性を効率的に上げるための研究を行っています。具体的には「多倍長数値計算(Multiple Precision Numerical Computation)」の研究が中心となります。現在のコンピュータが直接処理できるのは10進数で約16桁程度の有限桁小数ですが,どんな数も

$$ \frac{1}{3} \Longrightarrow 0.333333333333333 $$

と有限桁表現する必要があり,これを「丸め(rounding)」と呼びますが,この丸めによって生じる誤差が拡大することによって正しく計算されない「悪条件問題」というものがあり,シミュレーションの信頼性を損ねる一因となっています。それを防ぐ方法の一つとして,ソフトウェア的に桁数を増やして計算の精度を上げる,即ち,誤差を減らす手法が提案されており,それを「多倍長計算」と呼んでいます。

 多倍長計算の例として,単純なWebアプリケーションですが実行してみましょう。

Try MPFR!

 例えば,円周率\(\pi = 3.1415\cdots \)を10進約1000桁=2進約3321bitsで計算するには

try_mpfr_pi_3321

と入力すると

try_mpfr_pi_3321_result

という結果が得られます。答えが出力されているテキストボックス左下のところをドラッグすると拡大することができます。

 本研究室では,この多倍長計算を効率的に行うために,たくさんのパソコンをネットワークで繋いだPCクラスタを作ってきました。その際培われたネットワークのノウハウを生かし,Webアプリケーションを開発するようになり,現在は

  • 多倍長数値計算とその応用 = 高性能計算(HPC)
  • Webアプリケーションとの連携
  • 数学ソフトウェアの教育利用 → 中間的

にも勤しんでいます。ここで紹介する「MPFRgraph」はHPCとしての多倍長数値計算を広報する目的で構築されたものです。

3.Webアプリの仕組み

 ここで,Webアプリケーション(Web application),略して「Webアプリ」の構造を見ていくことにしましょう。先ほど紹介した「Try MPFR!」も一種のWebアプリです。これはもともと多倍長数値計算の土台を提供しているMPFRというライブラリ(沢山の機能を詰め込んだソフトウェアの素)を使って作った「mpfr_gexpr」という「ネイティブアプリ(native application)」をWeb経由で扱えるようにしたものです。

native_application_mpfr_gexpr

 ネイティブアプリは,パソコン(ハードウェア)を扱うための基本ソフトウェア(オペレーティングシステム, OS)の機能を利用して構築されたスタンドアローンタイプ(コンピュータ上で完結)のソフトウェアです。ハードウェア,OS,ネイティブアプリの関係性は,下記のようなレイヤー(層)で表現されます。

native_application

 これに対し,Webアプリは「Web (うぇぶ)」というインターネット上のサービスを土台としたアプリケーションで,「ブラウザ」と「サーバ」という,それぞれネットワーク上の独立したマシン上で協調して動作するものになっています。これらをレイヤー図で示すと下記のようになります。

web_application

 今のコンピュータは複数のプログラムを同時並行的に使用することができるので,ブラウザもサーバも同一マシンで動かすことは可能ですが,基本的にはインターネット上の独立したマシンでそれぞれ動くように作られています。ブラウザは今この記事をご覧になっている皆さんが直接目にしているソフトウェアです。代表的なブラウザとしては

があります。ちなみに,タブレットやスマートフォンではAndroid(Google)がChrome,iOSがSafariの利用が多いようです。Webサーバとしては

の利用が多いようです(参照: NetCraft Web server servey)。

 Webアプリの開発は,このようなブラウザとサーバの連携を前提にしてソフトウェア構築をする必要があります。そのためには次のようなコンピュータ言語(ソフトウェアを記述するための開発言語)が使用されます。

  1. HTML (HyperText Markup Language) ・・・ 静的Webページの構造を決める言語です。タグというもので括って文書構造を決定していきます。 → サンプル(注意:ビデオと音声が流れます)
  2. CSS (Cascading Style Sheet) ・・・ 静的Webページのデザイン(視覚効果)を決める言語です。
  3. JavaScript ・・・主としてブラウザ上で動作することの多いインタプリタ言語(一行単位で逐次解釈・実行される)です。
  4. PHP・・・主としてサーバ上で動作することの多いインタプリタ言語です。サーバ上のデータベースとの連携のためにも頻繁に利用されます。

 上記のうち,HTML, CSS, JavaScriptは一つの「HTMLファイル」(もしくは同等のもの)にまとめてサーバ側に保管され,24時間いつでもインターネット上のクライアントマシンのブラウザからアクセスできるようになっています。サーバ側はHTMLファイルの内容を送り出しているだけです。それがブラウザに届いて初めてレンダリング(画面構築)が行われ,皆さんが現在目にしているような画面として表示されているわけです(下記参照)。

web_application_process

 従って,ブラウザからユーザが入力したデータを受け取るには,サーバ側にもその仕組みが必要になります。ブラウザにデータ入力のためのインターフェースを作るためにはHTMLでフォーム(form)を記述することで可能になります。例としてアンケートフォームを示します。

アンケートフォーム(HTMLのみ)

 しかし,フォームを作っただけではサーバにデータは届きません。そのための受け皿として,PHPで記述されたプログラム,PHPスクリプトを作り,フォームと連携させる必要があります。上記のアンケートフォームにPHPスクリプトを連携させたものを下記に示します。

  1. アンケートフォーム(PHP連携:GETメソッド)
  2. アンケートフォーム(PHP連携:POSTメソッド)

 先に述べた,「Try MPFR!」もこれと同じ仕組みを使ってユーザからの情報「桁数(bits数)」と「数式(formula)」をサーバ側に伝えています。サーバ側ではそれをネイティブアプリである式パーサ「mpfr_gexpr」に与えて実行し,その結果をブラウザに送り込んでいます。

 この仕組みを拡張して関数の数表とグラフを同時に記述できるようにしたものが「MPFRgraph」です。

4.グラフ作成機能付き関数表作成アプリ「MPFRgraph」

 多倍長計算のためにはMPFRのような特殊なライブラリが必要ですので,Webアプリとしてネイティブアプリ多倍長精度式パーサmpfr_gexprを使って外部とのやりとりと連携させていました。しかしこの式パーサは数式に則って値の計算を行うことしかできません。多倍長計算のデモンストレーションとしては,様々な関数が使えることと,実際に関数 \(y = f(x) \) の右辺式を与えて,独立変数\(x\)を区間 \([\alpha, \beta]\)内で動かし,どのような値を取り,どのようなグラフになるのかを示すことが望ましいことになります。

 従って,区間\( [\alpha, \beta]\)と分割数\(n\),そして\(f(x)\)が\(x\)を含む数式として与えられていれば,分割された小区間幅を\( h = (\beta – \alpha) / n \)として与え,分点\(x_0, x_1, …, x_n\)を

\[ x_0 = \alpha,\ x_1 = x_0 + h,\ …,\ x_{n-1} = x_0 + (n-1) h,\ x_n = \beta \]

として与えることで,関数の値\(f(x_i)\) \( (i = 0, 1, …, n)\)とセットで

\[ (x_0, f(x_0)),  (x_1, f(x_1)), …, (x_n, f(x_n)) \]

をmpfr_gexprの機能を使うことだけで計算することができるようになります。こうして下記のように,グラフと数表を同時に出力できるようになりました。

 グラフについては,サーバ側でグラフィックスファイルととして表示させる方法と,ブラウザ側でHTML5のCanvas機能を使って表示させる方法があります。それぞれ作成してありますのでためしに様々な関数(多項式,三角関数,対数,指数等)を使って関数のグラフと数表を表示させてみてください。

 使い方はどれも同じです。例えば

\[ y = \sin x^2 + \exp(\cos x),\ x \in [0, 10],\ \mbox{分割数: }\ n = 10 \]

の場合は下記のように入力します。

input_example

 関数式を入力する際は四則演算,関数は下記のように表記します。これらを組み合わせ,変数名は必ず「x」として正しい数式を入力して下さい。

  • 加算:\( x + 3 \) → x + 3
  • 減算:\( x – 3 \) → x – 3
  • 乗算:\( 3x \) → 3 * x
  • 除算:\( \frac{x}{3} \) → x / 3
  • 平方根:\( \sqrt{x} \) → sqrt(x)
  • 三角関数: \( \sin x,\ \cos x,\ \tan x \) → sin(x), cos(x), tan(x)
  • 指数関数: \( e^x = \exp(x) \) → exp(x)
  • 自然対数: \( \log(x) \) → log(x)

 これらのアプリケーションの一連の処理の流れを示したものが下記の図になります。

 詳細は煩雑になりますので省略しますが

  1. 関数のグラフ描画するための関数値の計算は,CanvasバージョンではJavaScriptが担っている
  2. 多倍長精度の計算はサーバ側でしか実行できないので,関数表は全てサーバ側のmpfr_gexprが計算したものを使用している

という役割分担がなされています。

5.最後に

 高性能計算研究室ではほかにもいくつかのWebアプリや,Webプログラミング開発のための教材も公開しています。

幸谷研究室@静岡理工科大学

 また,「Webデザイン特別プログラム」では,過去の受講生の作品を公開しています。

Webデザイン特別プログラム

 24時間いつでも閲覧できますので,興味がありましたら色々見ていって下さい。特に優秀な作品は

Webデザイン特別プログラム優秀作品集

にまとめてあります。

参考URL

模擬講義「コンピュータはどうやって計算しているのか?」(改題「数学と計算 ~理論と実践~」)


概要

 数学(=理論)と計算(=実践)は全く性質の異なるもので,計算の効率性や正確性を追及すると数学とは別の方法論が必要になります。この講義では2次方程式や平均の計算を例に,「高性能な計算」をコンピューター上で行う手法を考えていきます。

1.数学=理論,計算=実践

 私の専門は「高性能計算(High Performance Computing)」です。読んで字のごとく,「高性能」な「計算法」を追求する学問の一分野ということになります。まだコンピューターが開発される前,20世紀までは人間があらゆる計算を行っていましたが,21世紀の現在,計算のうち単純なものはコンピューターが担うようになっています。従って,「高性能計算」も「コンピューターにとって」好都合なように計算方法を工夫するものになっています。この模擬講義では現在のコンピューターがどのような工夫のもとに「高性能計算」を実現しているかを具体的な例を挙げて解説します。

 さて,皆さんが思い浮かべる「計算」とはどのようなものでしょう? 「計算高い奴」と言うと,ずるがしこく立ち振る舞うという意味になりますが,ここで取り扱う計算とは「ごく単純な数の操作」とします。一番簡単な「数の操作」は,小学校で習う四則演算(加減乗除)ですね。ここで\(a\), \(b\)は実数とします。

  • たし算(加算)・・・\(a + b\)
  • ひき算(減算)・・・\(a – b\)
  • かけ算(乗算)・・・\(a \times b = ab\)
  • わり算(除算)・・・\(a / b\)

 四則演算だけで実行できる計算といえば,平均値の計算があります。例えば\(n\)個の実数\(c_1, c_2, …, c_n\)の平均値\(\mu\)は,\(n-1\)の加算と1回の除算を組み合わせて
\[ \mu = \frac{1}{n}\sum^n_{i=1} c_i = \frac{c_1 + c_2 + \cdots + c_n}{n} \]
と計算することができます。

 もう少し複雑な計算になると,平方根(\(\sqrt{a}\)),三角関数(\(\sin a, \cos a, \tan a \)),指数関数(\(\exp(x) = e^x\)),対数関数(\(\log a\))なども使う必要が出てきます。平方根を使う計算では,例えば2次方程式
\[ ax^2 + bx + c = 0 \]
の二つの解\(x_1\), \(x_2\)は
\[ x_1 = \frac{-b + \sqrt{b^2 – 4ac}}{2a}, x_2 = \frac{-b – \sqrt{b^2 – 4ac}}{2a} \]
として計算することができます。
 こういう計算は今では普通に電卓やパソコンソフト(電卓ソフトや表計算)で簡単に実行することができます。

 ところで,こういう計算はどうしてこういう式として表わすことができるのでしょうか? 「平均値」の計算はどうして全ての数を足してその個数で割るのでしょうか? 2次方程式の解はどうしてこのような形になるのでしょうか?

 そのような疑問に答えるための学問体系が「数学(Mathematics)」です。

 例えば平均値。解釈の仕方はいろいろ考えられますが,例えばこう考えてはどうでしょうか? いま,背丈の違う二人の人がいるとします。そのうち背の高い方から低い方へ背丈を分けて同じ身長にならしたとすれば,それがこの二人の「身長の平均値」を意味することになります。これが3人でも4人でも5人でも・・・\(n\)人でも,大きい人から小さい人へ背丈と分けて同じ身長に均した時の背丈が「\(n\)人の身長の平均値」であると言うことができます。

 2次方程式の解はどうでしょう? 教科書にはどうしてこういう解の公式が出てくるのか,ちゃんと解説が書いてあったはずですが皆さん覚えていますか?

 2次方程式の両辺を整理して,左側に\(x\)を含む1次式の2乗が来るようにすると
\[ \left(x + \frac{b}{2a}\right)^2 = \frac{b^2 – 4ac}{4a^2} \]
となりますので,
\[ x + \frac{b}{2a} = \pm \frac{\sqrt{b^2 – 4ac}}{2a} \]
から
\[ x = \frac{-b \pm \sqrt{b^2 – 4ac}}{2a} \]
となる訳です。

 こういう,「平均値の解釈」や「2次方程式の解の導出法」も「数学」の一部です。もちろん「平均値の計算」や「2次方程式の解の計算」も「数学」の一部ですが,これらの計算は「解釈」は「導出」の結果得られる「手続き」に過ぎません。計算の仕方が分かったとしても,どういう理由でそのような計算式になったのか,どういう理屈で計算式が導かれたのかは別に考える必要がある訳です。

 現在,コンピューターが担っている「計算」は数学的に導かれた解釈や導出法から得られる単純処理になっています。しかしこの先にもう一つ難しい問題が残っているのです。

 コンピューターはあらかじめ決められた手続きを教えてもらわなければ何も実行しません。手続きを記述した文書はプログラム(ソフトウェア)として人間が書いてやらなければなりません。平均値の計算や2次方程式の解の計算も,式に基づいてどういう手順で計算するのかを,四則演算や関数の組み合わせ方法まで指示しないと,効率的な計算ができないことがあります。このような手続きを「アルゴリズム(algorithm)」と呼びます。

 私が取り組んでいる「高性能計算」は,式を計算するためのアルゴリズムをどのように決定すれば効率が良くなるか,ということを考える学問分野なのです。

2.コンピューターの構造

 コンピュータは,全ての情報を電気信号として保持し,その信号を伝える回路(バス)と信号のやり取りを通じて何らかの処理を行う機器(デバイス)によって構築された金物(ハードウェア)と,ハードウェア内を伝わる電気信号の膨大な積み重ねによって構築された処理内容の蓄積,即ち軟体(ソフトウェア)によって構成されています。ハードウェアは堅い金属やプラスティックのカタマリですから,目で見ることができたり(可視物),触ってその存在を確認することができますが,ソフトウェアは単なる信号ですから,我々人間が感知するためにはデバイスを介して可視化する必要があります。処理を支持する信号を人間や外部の危機とやり取りするためのデバイスを入出力装置(I/O = Input & Output)と呼びます。こうした入出力装置を介しつつ,コンピューターに指令を与えるソフトウェアを作って初めてコンピューターを自在に操ることができるようになります。

basis_computing1

basis_computing2

 この原則は,先にあげた2次方程式の解の計算でも,平均値の計算でも変わりません。例えばC++(シープラスプラス)というコンピューター言語でこの計算を書くと次のようになります。これは計算部分のみの抜粋です(プログラム全体→quadratic_eq.cc)。

 解はx[0], x[1]という名前のメモリ上の領域に計算され格納されます。このプログラムはコンピュータが直接解釈できる形式に変換され(コンパイル),実行時に2次方程式の係数\(a\), \(b\), \(c\)を入力することで解\(x_1\), \(x_2\)のをディスプレイに表示(出力)します。例えば
\[ x^2 + 3x + 2 = 0 \]
の解\(-1, -2\)を計算する場合は次のように実行します。

$ ./quadratic_eq 1 3 2 ← プログラム名,係数の入力
1 * x^2 + (3) * x + (2) = 0 ← 方程式の確認
x[0] = (-1,0) ←解1
x[1] = (-2,0) ←解2

 また,平均を計算するC++プログラム”average.cc“は次のようになります。

 このプログラムをコンパイルして,\((1 + 2 + \cdots + 10) / 10 = 5.5 \)を計算させると,次のような結果が得られます。

$ ./average 10 1 2 3 4 5 6 7 8 9 10 ←プログラム名 データ数 データ1 ・・・ データ10
12
2
3
4
5
6
7
8
9
10
11
1 2 3 4 5 6 7 8 9 10 ←データの確認
average = 5.5 ←平均値

4.浮動小数点数と桁落ち
 さて,例として挙げた2次方程式の解計算と平均値の計算,実はどちらも正確なものではありません。コンピュータ内部での実数は,有限桁の小数でしか表現できないため,どうしても真の値とのズレ(誤差)が発生してしまいます。

 例えば,円周率\(\pi = 3.1415926535…\)は,同じ数のパターンの繰り返しが現れない無理数として良く知られたものですが,正確に表現するためには無限の桁数が必要になります。しかしコンピュータのメモリは有限なので,どこかで打ち切って丸めなければ格納できません。例えば15桁で表現するためには
\[ \pi \approx 3.14159265358979 \]
としなければならず,
\[\pi – 3.14159265358979 = 0.0000000000000032384626433832795… \]
という誤差が生じます。これを丸め誤差と呼びます。例として挙げた二つのプログラムではどちらも約15桁程度の正確さを持っておりますので,16桁目以降は不正確ということになります。

 実用的には15桁も正確であればさほど問題発生しませんが,稀に「桁落ち」という現象が起き,結果が著しく不正確になることがあります。例えば2次方程式の例では
\[ x^2 – 12345x + 0.00001 = 0\]
を計算してみると

$ ./quadratic_eq 1 12345 0.00001
1 * x^2 + (12345) * x + (1e-05) = 0
x[0] = (-8.10359779279679e-10,0)
x[1] = (-12344.9999999992,0)

という結果が得られますが,この場合\(x_1 = -8.100445524504379… \times 10^{-10}\)なので(正確な計算はこちらで実行できます),最初の3桁しか合っていません。誤差を含む,ほとんど同じ大きさの有限小数の加減算を行うとこのような「桁落ち」が起きます。これを防ぐには小数の桁数を増やすか,計算方法を変えるしかありません。

 ちなみに,桁数を約32桁に増やして計算してみると(プログラムはこちら


$ ./quadratic_eq_dd 1 12345 0.00001
1.000000000000000000000000000000e+00 * x^2 + (1.234500000000000000000000000000e+04) * x + (1.000000000000000000000000000000e-05) = 0
x[0] = (-8.100445524504379240908904588516e-10,0.000000000000000000000000000000e+00)
x[1] = (-1.234499999999918995544754956208e+04,0.000000000000000000000000000000e+00)

となり,桁落ち分をある程度フォローできていることが分かります。とはいえ,桁数を増やすより,もう少し賢いやり方はないのでしょうか?

 2次方程式の場合,次のように桁落ちしない正確な解\(x_1’\)を使って次の解\(x_2’\)を求めるという手順を取ることで,桁落ちを回避できることが知られています。
\begin{eqnarray*}
x_1′ &=& \frac{-b – \mbox{sign}(b)\sqrt{b^2 – 4ac}}{2a} \\
x_2′ &=& \frac{c}{2x_1′}
\end{eqnarray*}

 このような改良を加えたプログラムを”quadratic_eq2.cc“として実行してみます。

$ ./quadratic_eq2 1 12345 0.00001
1 * x^2 + (12345) * x + (1e-05) = 0
x[0] = (-12344.9999999992,0)
x[1] = (-8.10044552450438e-10,-0)

 このように,正確な解を得ていることが分かります。

 平均値の計算でも桁落ちが発生します。絶対値の小さい数から足しこむことである程度,桁落ちを防ぐことはできますが,2次方程式のほどの抜本的効果はありませんので,小数の桁数を増やすことが一番の解決策となります。

5.マルチコアと並列化

 コンピュータに求められる能力は増えることはあれ,減ることはありません。人間の利便性を高めるために高性能化は常に求められています。
 コンピュータの性能を高める方法は,CPUの処理能力を上げることです。その一つの方法として,動作周波数(clock)を高めるという方法があります。実際に,20世紀最後の10年ほどで飛躍的に動作周波数は増加しました。

Intel CPUの動作周波数

 しかし,21世紀に入ると打ち止め状態になっています。動作周波数を高めると電力を食うようになります。一般家庭でも使用できるようなパソコンが何千ワットも食うようでは使い物になりません。また,電力消費に伴う発熱も深刻で,安定的に動作させるための空冷ファンは大きくなるばかり。余計に電力消費が増えていきます。
 性能は上げなくてはならない,しかし,電力消費量の増加は極力抑えたい・・・という二律背反な要求に対し,CPUメーカーは電子回路の改良と共に,「マルチコア」という方向にかじを切りました。簡単に言うと,CPUの処理の中核を担う回路(コア)を複数(マルチ)持たせて,あたかも複数のパソコンを同時に使えるようにする技術です。もし6つのコアを持つCPUがあれば,6つの処理を同時にこなすことができ,短時間で終了することで多少電力消費は増えてもトータルでは少なくすることができるはずです。

muticore

 このように,複数の処理を同時並行にこなすことを「並列処理」と呼びます。対して,一つの処理だけを順次行う処理を「逐次処理」と呼びます。現在のパソコンやスマートフォン(スマホ)は2~4コアCPUが普通ですし,GPUのようにもっと大量のコアを備えてグラフィックスや計算を高速にこなす「メニー(many, 多数)コア」ハードウェアも登場しています。

gpu

 では,どうやったら複数のコアを当時に使うことができるのでしょうか? 並列処理を行うプログラム技法はいくつかありますが,一つのプログラムの中の処理を「スレッド」という単位に分割して並列に実行可能にさせる,というマルチスレッド化が良く用いられます。例えば,OpenMPという技法を使って2次方程式の解を同時に計算させるプログラムは次のようになります。

\(x_1, x_2\)の計算は相互に依存していませんので,並行して計算させることができるわけです。逆に,桁落ちが起きないように改良してしまうと,片方の解が出るまでは次の解の計算に進めず,このような並列処理はできなくなります。

 しかし,2次方程式の計算は大変軽いものなので,並列化しても大したメリットはありません。桁落ちを防ぐ処理を行う方が実用的には有用と言えるでしょう。

6.並列処理の効果のほどは?

 では平均の計算はどうでしょう? 大量のデータ(数万,数億,・・・)の平均値の計算は普通に行われるものなので,並列処理するメリットはかなりありそうです。ではどのように行えばいいのでしょうか? 前述した

\[ \frac{1 + 2 + \cdots + 9 + 10}{10} = 5.5 \]

の計算を例に考えてみましょう。

 2コア(デュアルコア)のCPUでは,同時に2スレッドの処理が可能です。その際には,半分ずつ和の計算を行い,次にその部分和を足し合わせて1スレッドで合計値を求め,最後に10で割ります。

parallel2

 4コア(クアッドコア)が使える場合はどうでしょう? 最初の部分和は4スレッドで計算し,次の部分和の計算は2スレッドで分担します。何故なら,最初の部分和の計算よりも計算量が増えると,ここで時間が余計にかかってしまうからです。

 最後は2個の和の計算を1スレッドで行って10で割ります。

parallel4

 この場合,和の計算は最大10回で済みますから,5スレッド以上使える環境があっても意味がありません。

 こうして2スレッド,4スレッドそれぞれ並列実行できたとすると,計算時間は4スレッドの方が短くなることが期待されます。

parallel_time

 では更に大量のデータの平均の計算を実行しようとするとどうでしょうか? データ数\(n\)に対して,同時に実行可能なスレッド数を\(\alpha\) (\(n > \alpha\))とすると最初の部分和の計算は\(\lceil n / \alpha\rceil \)回の和で済みます。ここで\(\lceil x \rceil\)は \(x\)を超えない最大の整数を意味します。

 次の段階で足すべきデータは\( n/2 \)個になっています。これを繰り返すことで,最後は2個の和で済ますことができます。最終的に2個になったデータの和と\(n\)の除算は1スレッドで終了しますので,トータルでは和の計算の部分は最大
\[ \lceil \log_2 n \rceil \]
回の繰り返しで平均の計算が完了することになります。先にあげた10個のデータの場合は\(\lceil\log_2 10\rceil\) \(=\lceil 3.321928…\rceil = 4\)となり,和の計算は4段階必要なことになります。

 このような集約を行う計算をリダクションと呼びます。これをOpenMPで記述すると下記のようになります。(プログラム全体: average_omp.cc)

 使用するスレッド数はomp_set_num_threads関数で設定しています。データの数が少ない時はスレッド数を制限するようにしています。

 では実際に並列処理によってどの程度計算時間が軽減されるかを見ていきましょう。パソコンの環境は次の通りです。

  • CPU: Intel Core i7-3820 3.6GHz 4コア
  • メモリ: 32GB
  • OS: CentOS 6.5 x86_64
  • コンパイラ: Intel C++ 14.0.2
  • データ量: 16GB (\(= 8\mbox{Bytes} \times 2 \times 1024^3\))

4コアありますので,最大4スレッドを並列に実行できることになります。そこで,スレッド数を1, 2, 4として実行してみました。”Time(s)”の所が処理時間(秒)になります。


$ export OMP_NUM_THREADS=4;./average_omp 1
#threads: 1
num_c : 2147483648
average = 1073738907.87731
Time(s) = 0.998027801513672
$ export OMP_NUM_THREADS=4;./average_omp 2
#threads: 2
num_c : 2147483648
average = 1073738907.87725
Time(s) = 0.568727016448975
$ export OMP_NUM_THREADS=4;./average_omp 4
#threads: 4
num_c : 2147483648
average = 1073738907.87729
Time(s) = 0.453949928283691

1スレッドでは約1秒でしたが,2スレッドで約0.57秒,4スレッドで0.45秒で終了しています。2スレッドで大体半減しますが,4スレッドにすると更に半分に,とはなりません。スレッド数が多いとその準備の手間で,理論的な見積りよりも処理時間が余計にかかります。

 ちなみに,太線部分のところがスレッド数によって異なっていることにお気づきでしょうか? 並列化してしまうと,特段指定しない限り,計算の順番は全く保障されないので,丸め誤差の差異が発生しているのです。これは桁落ちしていないのでせいぜい末尾の2桁程度で済んでます。

7.まとめ

 以上,駆け足でしたが,数学と計算の違い,そして計算そのものもきちんと考えることは意外と面倒なことが理解できたことでしょう。今まで述べてきたことのポイントを下記にまとめておきます。

     

  • 理論式をそのまま逐次処理的に計算するにしても,桁落ちが起きないような工夫が可能であれば行う必要があります。
  • マルチコアCPUでは並列処理が可能なので,大量のデータの処理を行うときには威力を発揮できます。しかしそれは並列処理可能な部分に限られます。
  • 並列処理を行って性能向上を図るためには,処理の流れ,即ちアルゴリズムを工夫する必要があります。

 現在の科学技術計算は,シミュレーションの理論化,計算処理,それぞれ分担して専門家が担うようになっています。ここで述べたのはごく簡単な例に基づいた計算処理の部分だけですが,数式をそのまま計算するというだけでは済まない面倒さがあることを覚えて頂ければ幸いです。

京都大学のOffice365サービス移行記録

 本学は本年度(2016年度)よりMicrosoftのOffice365サービスの利用を開始しました。メールが世界どこにいても読み書き返信できるようになり,Office 2016の利用も手軽になるなど良いことづくめのようですが,反面,すべての情報を先方に預けるということでもあり,良いことづくめというわけではありません。

 京都大学は本学より早くOffice365とその前身となるサービスへの移行を開始しました。その貴重な記録が公開されています。ネットワーク管理やOffice365サービスに興味のある方は読む価値があります。

 ネットワーク管理者としては,サーバのメンテナンス業務が楽になる反面,こちらではどうしようもないトラブルへの対処も先方任せになるというジレンマに見舞われることになる訳で,日々胃が痛いことには変わりないなぁという感想を持ちました。

「効率的な手順とは?」 : 2/1 裾野高校模擬講義

1.初めに
 世の中にはいろいろな「作業」があります。楽しい作業,苦しい作業,大変な作業,楽な作業・・・傍から見ているのと実際にその作業を行ってみるのでは全く異なる感想を持つ,なんてことはよくあることです。「計算」なんてその最たるものでしょう。数学と聞くともう数式を見るのも嫌い!,という人は多いですが,小学校の時に習った四則演算(足し算,引き算,掛算,割算)であれば,楽しかったかどうかはともかく,少ない桁数の自然数どうしでしたらそれほど時間をかけることなく実行できるはずです。例えば\[2\times 2\times 2\times 2\times 2 = 2^5 = 32\]なんて計算ができないという人はまずいないでしょう。

 さて「作業」にもイロイロありますが,ここではその開始から終了までに要する「時間」を考えてみたいと思います。作業に要する時間は,どのぐらいの手間がかかるのか,つまり,その作業を行うために必要となる最も簡単かつ短時間で終わる小作業の回数をカウントすることで大体の見積もりを取ることができます。先に挙げた\(2^5\)を計算する作業は,2を掛けるという小作業を4回繰り返すことで32という結果を得て完了することができました。従って,これが\(2^6\), \(2^7\), …, \(2^n\)であっても,掛算をそれぞれ5回,6回,…, \(n\)回実行すれば,計算結果を得ることができ,作業を完了することができるわけです。

 しかしもっとこの作業を効率的にできないものでしょうか? 同じ作業を実行するにしても,初めて経験する人よりは,普通は,その作業を経験したことのある人の方が「効率的」=「短い時間で済む」=「少ない小作業の回数で済ませる」ことができますよね? それは同じ作業を繰り返すことで,「この小作業は不要,ここでこの小作業をすることで効率的になる」ということに気が付くからです。現在はどんな工場でも,どんな事務所でも,どんな会社でも,短い時間で効率的に作業を行うことが求められます。これができないとどうなるか? そう,同じ仕事量でも完了までにいたずらに時間がかかり,いつまでたっても帰れない「ブラック企業」と言われることになってしまうわけです。サービス残業させられて残業代も出ないようなら,短い時間で済ませてさっさと定時退社した方がいいですよね?

 さて,先ほど実行した\(2^n\)の計算ですが,もっと効率的に実行する方法がありますよね? そう,掛算の回数を減らすことはそんなに難しくないですよね? まぁ電卓が手元にあれば,という前提で考えてもらうと,例えば\[2^5 = (2\times 2)^2 \times 2 = 32\]という手順であれば,3回の掛算で済むことになります。\(2^6\)ですと,\[ (2\times 2)^3 = (4\times 4) \times 4 = 64\] となってこれも3回で済みます。・・・と考えていくと,一般に\(2^n\)の場合,\(n\)の2進表記が\(n = (p_k p_{k-1} \cdots p_1 p_0)_2\)であれば,\[2^n = 2^{2^k\cdot p_k} \times 2^{2^{k-1}p_{k-1}} \times \cdots \times 2^{2{p_1}}\times 2^{p_0} \]という手順で済むことが分かります。つまり\(2^{100}\)であれば,99回の掛算が必要になるところを,\(100 = (110010)_2\)より,\[ 2^{100} = 2^{2^6+2^5+2^2} = 2^{64 + 32 + 4} = 2^{64}\times 2^{32} \times 2^4 \]となり,\[ 2^{64} = (((((2^2)^2)^2)^2)^2)^2 \] ですから6回の掛算で済み,この途中結果を覚えておくとすると自動的に\(2^{32}\)と\(2^4\)が得られますので,合計8回の掛算で済むわけです。ちなみに答えは\[\begin{split} 2^{100} &= 18446744073709551616\times 4294967296\times 16\\ &= 1267650600228229401496703205376\end{split} \] となります。

 単純な掛算でも計算回数を減らすことができるわけです。これが膨大な数の掛算を行うとなると,効率的な「手順」=「アルゴリズム」を使わない手はありません。

2.ハノイの塔問題
 ではもう少しゲーム的な「作業」を行ってみましょう。アルゴリズムを考えるときには必ず題材になる有名なパズル「ハノイの塔」というものです。
 次のような積木があるとします。

hanoi_3units

 この積木を一番端の棒に移動させて下さい。但し次の条件を必ず守って下さい。

  • 円盤は1枚ずつしか移動できません。移動先は必ず他の棒になります。
  • 小さい円盤の上に大きな円盤を載せてはいけません。
  • 作業は全ての円盤が一番端の棒に,大きさ順に刺さったところで完了します。

 さて,ではこの3枚のハノイの塔問題を,なるべく「効率的」に解いて下さい。

 解答は下記のYouTube動画にあります。

3.トランプの並べ替え:クイックソート
 日常的に良くある作業として「並べ替え」があります。これを効率的にできないものでしょうか?

 例えばジョーカーの入っていない52枚のトランプをスペード,クローバー,ハート,ダイヤモンドの組に分け,それぞれ番号の小さい順に並べ替える問題を考えてみます。どうすればなるべく早くに並べ替えが完了するでしょうか?

 例えば,下記のYouTube動画に示す「クイックソート」という手順は如何でしょうか?

4.おわりに
 コンピューターは人間と違って莫大なデータ量(ビッグデータ)を扱うために使用されます。従ってちょっとでも非効率的なアルゴリズムを使うと途端に膨大な時間を無駄にすることになりかねません。効率的な手順を考えることは,コンピューターを効率的に使うことにもつながります。作業が早く完了すればそれだけ電気代が節約できますからね。
 我々コンピューターの専門家はこのような「効率的」なアルゴリズムを考えることが重要な仕事の一つなのです。

[模擬講義2] 大規模計算とCPU&GPU


 これは模擬講義「数学と計算 ~理論と実践~」の続きになります。まずそちらを読んで頂いてから,こちらの文章に取り組んで下さい。

1. 大規模計算とは?

 大規模計算とは,現在のコンピュータが実行するに際し,大量のデータをメモリに格納し,計算時間にも数秒,数分,数時間・・・数日・・・費やす必要のある一連の計算処理のことです。
 既にお話しした通り,コンピュータは,処理する対象となるデータを一度メモリに格納し,CPUでそれらを読み込んで処理し,処理結果をメモリに書き戻すという作業を繰り返します。この作業1回分の処理時間はごく短く,マイクロ秒(μs)\(=10^{-6}\)秒単位の世界です。しかし例えば5μsの処理を10回行えば50μs,100回行えば500μs,…と,繰り返しの数\(n\)が大きくなればなるほど,合計処理時間\(5n\) μsは大きくなります。このように\(n\)に比例して時間が増える処理のことを「\(O(n)\)のアルゴリズム」と称します。先に挙げた平均の計算
\[ \mu = \frac{1}{n} \sum^n_{i=1} c_i = \frac{c_1 + c_2 + \cdots + c_n}{n} \]
は典型的な\(O(n)\)のアルゴリズムです。

 しかし,実用的な処理の場合,規模,即ち\(n\)が大きくなるにつれてもっと大きな処理時間を要するものが存在します。少し話が難しくなりますが,以下では連立一次方程式(連立方程式の一種)を見ていくことにしましょう。

2. 連立一次方程式と大規模化

 実数の定数\(a\), \(b\), \(c\), \(d\)に対して,次の二つの方程式を同時に満足する解\(x\), \(y\)を求める,というのが,中学校で初めて習う連立方程式,正式には連立「一次」方程式です。

\[\left\{\begin{array}{ccc} ax + by &=& \alpha \\
cx + dy &=& \beta \end{array}\right. \]

 具体的には例えば

\[\left\{\begin{array}{ccc} 3x + 2y &=& 1 \\
-4x – 2y &=& -2 \end{array}\right. \]

のようになります。この場合,\(x = 1\), \(y = -1\)が解となります。

 連立一次方程式は,様々な場面で登場するとても重要な方程式です。そしてそれだけに,とても大規模なものを扱うことが良く行われています。そのため,大規模化しやすいよう,定数や変数は\(a\)と\(x\)という二つの文字に制限し,大規模化した時には定数や変数に添え字(下付き数字)をつけて区別するように変更します。例えば上の2元2次連立一次方程式は下記のように書き直します。

\[\left\{\begin{array}{ccc} a_{11}x_1 + a_{12}x_2 &=& b_1 \\
a_{21}x_1 + a_{22}x_2 &=& b_2 \end{array}\right. \]

 添え字が異なる文字は違うものであることを意味します。ちょっと目がチカチカしますが,我慢して下さい。

 さて,連立一次方程式の重要なところは,変数の個数と方程式の個数が同じであるところです。2元2次の場合は,それぞれの方程式が直線になっており,2つの直線の交わる交点が解となっている,ということを考えると,解が一つに定まるためには2つの方程式=2本の直線が必要であることは感覚的に理解できるでしょう。

 では変数が3つ,即ち\(x_1\), \(x_2\), \(x_3\)になった時はどうでしょう? この場合,方程式は
\[ ax_1 + bx_2 + cx_3 = d \]
という形になりますので,図形としては平面を意味することになります。つまり「板(いた)」ですね。並行でない「板1」と「板2」が二つ重なるところは直線になりますので,点としては一つに絞れません。更にもう一つ「板3」を持ってくると,既に存在している「板1」と「板2」が交わる「直線12」がこの「板3」とぶつかる交点ができますので,これが解として一つに定まることになります。つまり,一つの解=交点が決まるには,板=方程式は3つ必要になります。従って3つの変数に対しては下記のような3元3次連立一次方程式である必要がある訳です。

\[\left\{\begin{array}{ccc} a_{11}x_1 + a_{12}x_2 + a_{13}x_3 &=& b_1 \\
a_{21}x_1 + a_{22}x_2 + a_{23}x_3 &=& b_2 \\
a_{31}x_1 + a_{32}x_2 + a_{33}x_3 &=& b_3 \\
\end{array}\right. \]

 以下は同様です。変数が4つになれば方程式は4本必要になります。変数が5つになれば方程式は5本,・・・,つまり,変数が\(n\)個あれば方程式は\(n\)本必要になる訳です。それを下記のように表現します。

\[\left\{\begin{array}{cccc} a_{11}x_1 + a_{12}x_2 + \cdots + a_{1n}x_n &=& b_1 \\
a_{21}x_1 + a_{22}x_2 + \cdots + a_{2n}x_n &=& b_2 \\
\cdots & & \cdots \\
a_{n1}x_1 + a_{n2}x_2 + \cdots + a_{nn}x_n &=& b_n \\
\end{array}\right. \]

 さらに簡略化して書きたいときには,係数\(a_{ij}\),変数(解)\(x_i\),右辺の定数項\(b_i\)をバラシて
それぞれのカタマリを,\(A\), \(\mathbf{x}\),\(\mathbf{b}\)と表記し

\[ A = \left[\begin{array}{cccc}
a_{11} & a_{12} & \cdots & a_{1n} \\
a_{21} & a_{22} & \cdots & a_{2n} \\
\vdots & \vdots & & \vdots \\
a_{n1} & a_{n2} & \cdots & a_{nn}
\end{array}\right],\ \mathbf{x} = \left[\begin{array}{c}
x_1 \\
x_2 \\
\vdots \\
x_n
\end{array}\right],\ \mathbf{b} = \left[\begin{array}{c}
b_1 \\
b_2 \\
\vdots \\
b_n
\end{array}\right] \]

以下のように一行の数式で表現します。

\[ A\mathbf{x} = \mathbf{b} \]

ここで\(A\)を係数からなる「行列(matrix)」,\(\mathbf{x}\)を「解ベクトル(vector)」,\(\mathbf{b}\)を「定数ベクトル」と呼びます。最初に例示した2元2次方程式

\[\left\{\begin{array}{ccc} 3x_1 + 2x_2 &=& 1 \\
-4x_1 – 2x_2 &=& -2 \end{array}\right. \]

に対しては
\[ A = \left[\begin{array}{cc} 3 & 2 \\
-4 & -2 \end{array}\right],\ \mathbf{x} = \left[\begin{array}{c}
x_1 \\
x_2
\end{array}\right],\ \mathbf{b} = \left[\begin{array}{c}
1 \\
-2
\end{array}\right] \]
となります。

3.連立一次方程式の解き方

 さて,では大規模な連立一次方程式をどう解いていくか? その方法を説明していきます。詳細については大学学部の「線形代数」という科目で習う内容になるので,そのエッセンスだけを述べることにします。

 まず,簡単な2元2次の例で確認してみましょう。肝心なのは,「自動的に解ける形に方程式を変形する」→「自動的に解いてしまう」という手順を踏むということです。例えば\(u_{21} = 0\)であれば

 \[\left\{\begin{array}{rcl} u_{11}x_1 + u_{12} x_2 &=& y_1 \\
u_{22} x_2 &=& y_2 \end{array}\right. \]

と表現することができますので,下の式から順に\(x_2 = y_2 / u_{22} \)として\(x_2\)を求め,次にそれを上の式に代入して\(x_1 = u_{12} x_2 / u_{11}\)を求めればよいことになります。このように,一番下の式が変数が一つだけの方程式になった格好の連立一次方程式は,係数行列\(A\)が上三角行列\(U\),即ち
\[ U = \left[\begin{array}{cc}
u_{11} & u_{12} \\
0 & u_{22}
\end{array}\right] \]
を用いて
\begin{equation}
U\mathbf{x} = \mathbf{y} \label{eqn:upper_triangle}
\end{equation}
と表現することができます。このように上三角行列を係数として持つ連立一次方程式を下の方から解いていく方法を「後退代入」と呼びます。

 同様に\(l_{11} = l_{22} = 1\)かつ\(l_{12} = 0\)であれば

\[\left\{\begin{array}{rcl} y_1 \phantom{+ u_{12}y_2} &=& b_1 \\
l_{21}y_1 + y_2 &=& b_2 \end{array}\right. \]

と表現できますので,これも上から順に\(y_1 = b_1\),次に\(y_2 = b_2 – l_{21}y_1\)として解を求めることができます。このとき係数行列は下三角行列\(L\),即ち
\[ L = \left[\begin{array}{cc}
1 & 0 \\
l_{21} & 1
\end{array}\right] \]
を用いて
\begin{equation}
L \mathbf{y} = \mathbf{b}
\end{equation}
と表現することができます。このような下三角行列を係数として持つ連立一次方程式を上の方から解いていく方法を「前進代入」と呼びます。

 結論を先に述べると,上記2つの形の連立一次方程式を導いて,次のようにして解きます。

  1. [LU分解]  係数行列\(A\)から下三角行列\(L\)と上三角行列\(U\)を求める(実は\(A = LU\)という関係になります)
  2. [前進代入] \(L\mathbf{y} = \mathbf{b}\)を\(\mathbf{y}\)について解く
  3. [後退代入] \(U\mathbf{x} = \mathbf{y}\)を\(\mathbf{x}\)について解く

 実際に
\[\left\{\begin{array}{ccc} 3x_1 + 2x_2 &=& 1 \\
-4x_1 – 2x_2 &=& -2 \end{array}\right. \]
でLU分解→前進代入→後退代入を行っていきます。

 まずLU分解です。
\[ A = \left[\begin{array}{cc}
3 & 2 \\
-4 & -2
\end{array}\right] \]
を,行(横一列)をひとまとめに足したり引いたり定数倍したりして,上三角行列\(U\)にします。
\[ \left[\begin{array}{cc}
3 & 2 \\
-4 & -2
\end{array}\right] \Rightarrow \left[\begin{array}{cc}
3 & 2 \\
-4 – (-4)/3\times 3 & -2 – (-4)/3\times 2
\end{array}\right] = \left[\begin{array}{cc}
3 & 2\\
0 & 2/3 \\
\end{array}\right] = U \]
第1行目に\((-4)/3\)を掛けて,第2行目から引いているわけです。こうして\(a_{21}\Rightarrow 0\)にすることができました。実はこの時,下三角行列\(L\)も自動的に求められており,
\[ L = \left[\begin{array}{cc}
1 & 0 \\
-4/3 & 1
\end{array}\right] \]
になります。実は\(a_{21} / a_{11} = -4 / 3\)を計算しただけになります。

 こうして出てきた\(L, U\)を用いて
\[ L\mathbf{y} = \mathbf{b} \Longleftrightarrow \left[\begin{array}{cc}
1 & 0 \\
-4/3 & 1
\end{array}\right] \left[\begin{array}{c}
y_1 \\
y_2
\end{array}\right] = \left[\begin{array}{c}
1 \\
-2
\end{array}\right] \]
を解いて\[ \mathbf{y} = \left[\begin{array}{c}
1 \\
-2/3
\end{array}\right] \] を求め,次に
\[ U\mathbf{x} = \mathbf{y} \Longleftrightarrow \left[\begin{array}{cc}
3 & 2 \\
0 & 2/3
\end{array}\right] \left[\begin{array}{c}
x_1 \\
x_2
\end{array}\right] = \left[\begin{array}{c}
1 \\
-2/3
\end{array}\right] \]
を解いて\[ \mathbf{x} = \left[\begin{array}{c}
1 \\
-1
\end{array}\right] \] を求めます。

 この方法は\(n\)元\(n\)次連立一次方程式になっても変わりません。しかし,計算の手間が膨大に増えていきます。実際,LU分解は
\[ A = \left[\begin{array}{cccc}
a_{11} & a_{12} & \cdots & a_{1n} \\
a_{21} & a_{22} & \cdots & a_{2n} \\
\vdots & \vdots & & \vdots \\
a_{n1} & a_{n2} & \cdots & a_{nn}
\end{array}\right]\]
\[\Rightarrow \left[\begin{array}{cccc}
a_{11} & a_{12} & \cdots & a_{1n} \\
a_{21}/a_{11} & a_{22} – a_{21}/a_{11}\times a_{12} & \cdots & a_{2n} – a_{21}/a_{11}\times a_{1n} \\
\vdots & \vdots & & \vdots \\
a_{n1}/a_{11} & a_{n2} – a_{n1}/a_{11}\times a_{12} & \cdots & a_{nn} – a_{n1}/a_{11}\times a_{1n}
\end{array}\right]\]
\[\Rightarrow \cdots \Rightarrow \left[\begin{array}{cccc}
a_{11} & a_{12} & \cdots & a_{1n} \\
a^{(1)}_{21} & a^{(2)}_{22} & \cdots & a^{(2)}_{2n} \\
\vdots & \vdots & & \vdots \\
a^{(1)}_{n1} & a^{(2)}_{n2} & \cdots & a^{(n-1)}_{nn}
\end{array}\right]\]
となります。この過程を図示すると下記のようになります。

lu

黒く塗りつぶした四角部分は
\[ a^{(i-1)}_{ji} / a^{(i-1)}_{ii}\ (j = i + 1, …, n) \]
の計算を行い,灰色の部分は
\[ a^{(i-1)}_{jk} – a^{(i-1)}_{ji} / a^{(i-1)}_{ii} \times a^{(i-1)}_{ik}\ (j, k = i+1, …, n) \]
の計算を行っています。

黒い四角部分は\[ (n-1) + (n-2) + \cdots + 2 + 1 \]個となり,灰色の部分は\[ (n-1)^2 + (n-2)^2 + \cdots + 2 + 1 \]個となり,それぞれ1回ずつの割算と掛算が入っていますので,この割算と掛算の回数だけでも
\[ \sum^{n-1}_{i=1} i(i + 1) = \frac{n(n^2 – 1)}{3} \]
となりますから,最終的には\(O(n^3)\)の手間がかかるようになります。\(n\)が10倍ずつ増えていくとき,\(O(n)\)の計算が10秒,100秒,1000秒・・・と増えていくとすると,\(O(n^3)\)は1000秒,1000000秒, 1000000000秒・・・のように増えていきます。計算時間の増え方が1000倍になっていくわけです。

 こうして\(L, U\)は
\[ L = \left[\begin{array}{cccc}
1 & 0 & \cdots & 0 \\
a^{(1)}_{21} & 1 & \cdots & \vdots \\
\vdots & \ddots & \ddots & 0 \\
a^{(1)}_{n1} & \cdots & a^{(n-1)}_{n,n-1} & 1
\end{array}\right],\ U = \left[\begin{array}{cccc}
a_{11} & a_{12} & \cdots & a_{1n} \\
0 & a^{(2)}_{22} & \cdots & a^{(2)}_{2n} \\
\vdots & \ddots & & \vdots \\
0 & \cdots & 0 & a^{(n-1)}_{nn}
\end{array}\right] \]
として得られますので,前進代入は\(L\mathbf{y} = \mathbf{b}\)を解いて
\begin{align*}
y_1 &:= b_1 \\
y_2 &:= b_2 – a^{(1)}_{21} y_1 \\
& \vdots \\
y_{n-1} &:= b_{n-1} – \sum^{n-2}_{j=1} a^{(j)}_{n-1,j}y_j \\
y_n &:= b_n – \sum^{n-1}_{j=1} a^{(j)}_{nj}y_j
\end{align*}
となり,後退代入は\(U \mathbf{x} = \mathbf{y}\)を解いて
\begin{align*}
x_n &:= \frac{y_n}{a^{(n-1)}_{nn}} \\
x_{n-1} &:= \frac{y_{n-1} – a^{(n-2)}_{n-1,n} x_n}{a^{(n-2)}_{n-1, n-1}} \\
& \vdots \\
x_2 &:= \frac{y_2 – \sum^{n}_{j=3} a^{(1)}_{2j} x_j}{a^{(1)}_{22}} \\
x_1 &:= \frac{y_1 – \sum^{n}_{j=2} a_{1j} x_j}{a_{11}}
\end{align*}
として計算することができます。

 従って,前進代入,後退代入の手間は\(O(n^2)\)になります。

4. 並列で計算するには?

 さて,このように大規模化した連立一次方程式を現在のコンピュータ上で高速に解くためには,LU分解,前進代入,後退代入の計算を並列化する必要があります。

 LU分解は下記のように並列化します。色分けした部分を各スレッドに分配して並列に計算しています。

parallel_lu

 前進代入,後退代入も下記のように並列化します。

parallel_subst

 LU分解も前進代入,後退代入も,計算が進むにつれて計算量が減り,結果的にだんだん並列化効率が落ちていきますが,これ以上の良い手段はなかなか見つからないというのが現状です。しかしこのように並列化することで,マルチコアCPU環境・GPU環境では高速に計算することができるようになります。実際に12コアCPUとTesla K20というGPUを積んだワークステーション環境で解いた結果が下記になります。

dsgesv_cpu_gpu

 最大\(n = 8192\)までの規模の連立一次方程式を説いた時の計算時間と計算に伴う丸め誤差をプロットしています。GPUの方が圧倒的に高速に計算できていることが分かります。

5. まとめ

 駆け足で解説してきましたが,ここで説明した大規模な連立一次方程式をLU分解・前進代入・後退代入を使って解く方法は直接法(direct method)と呼ばれ,Linpack(リンパック)ベンチマークテストと呼ばれる,スーパーコンピュータの速度比較プログラムに使われています。Top500というサイトには,毎年半期ごとに世界中のスーパーコンピュータのLinpackテストの結果に基づいたランキングが掲載されており,2015年9月現在のナンバーワンは中国にあります。日本では神戸にある「京(ケイ)」というスーパーコンピュータが最高速で,現在は世界第4位に後退しました。

 コンピューターの性能評価のためのLinpackテストは,皆さんが中学校で習う連立一次方程式の拡大版で行われており,それらは\(O(n^3)\)の処理時間を要するものである,ということを理解してもらえば幸いです。