パズルを解く制約プログラミング

このような10枚の板を、四角に組み合わせるというパズル。

行き当たりばったりでは、簡単には解けない。
パソコンに解かせようと探したところ、こんなプログラムを見つけた。

* Coprisによる制約プログラミング入門
>> http://bach.istc.kobe-u.ac.jp/copris/docs/intro-ja.html

“制約プログラミング”とは、制約条件をコンピュータに入れると、コンピュータが条件に合った答をはじき出す、というもの。
Coprisの場合、条件を整数の数式の形で入力すると、答にあてはまるパターンが次々と出てくる。
この感覚を口で説明するのは、なかなか難しい。
たとえば上のパズルを解くプログラムは、こんな感じになる。

* まず答が入る容れ物となる、5x5の変数を用意する。
 この全部で25個の変数は、それぞれ -1, 0, +1 の値のいずれかを取るものとする。
 +1 は長い切れ込み、-1 は短い切れ込み、0 は普通の長さの切れ込みを表すことにする。

* 10枚の板には +1 と -1 が1個ずつ入っている。
 [条件1] 5x5の変数の、それぞれ縦の列の合計 = 0 となる。
 [条件2] 5x5の変数の、それぞれ横の行の合計 = 0 となる。
 [条件3] 5x5の変数の、それぞれ縦の列の絶対値の合計 = 2 となる(+1 と -1 の2個が入っている)。
 [条件4] 5x5の変数の、それぞれ横の行の絶対値の合計 = 2 となる(+1 と -1 の2個が入っている)。

* 10枚の板のパターンが全部異なっている。
 板を裏返しに差し込むこともできるが、それらのパターンも全て異なる。
 [条件5] 5カ所の切れ込みを2進数と見なしたとき、2進数の値が縦横全部で異なる。
  2進数は、切れ込みを右から左へ読んだパターンと、左から右に読んだパターンの全てが異なる。

* Scalaソースコードはこちら・・・
>> http://brownian.motion.ne.jp/memo/Copris/CrossBoard.scala
>> http://brownian.motion.ne.jp/memo/Copris/CrossBoardMain.scala

このような制約条件をセットして、答を探せ(find)と命令すると、Coprisが次々と答を出してくる。
試したところ、176パターンの答が出てきた。
このパターンの中には上下左右前後をひっくり返しただけの答も含まれているので、実質的には 176÷8=22 パターンの答があるようだ。
1パターンだけ示すと、こんな風になる。

0 -1 0 0 1
-1 1 0 0 0
0 0 1 -1 0
1 0 0 0 -1
0 0 -1 1 0

* 全ての解はこちら・・・
>> http://brownian.motion.ne.jp/memo/Copris/CrossBoardSol.txt

■ 班分け問題

パズルを解くプログラムなんて、ピンポイントでマニアックなものかと思いきや、これが案外役に立つ。
実際、私の役に立ったのは“相性のあるグループ分け問題”だった。
40名ほどのメンバーを、5〜6名×7グループに班分けしたかったのだが、
メンバー同士には相性があって、この人と組みたい、この人とは一緒になりたくない、といった希望がある。
これが40人分ともなると、いちいち希望を聞いて班分けするのは実に面倒くさい。
場合によっては、あっちの希望はかなったのに、なぜこっちの希望はかななわないのか、など、不平不満になりかねない。
そこで制約プログラミングの出番である。

* まず答が入る容れ物となる、(メンバー数)x(グループ数)の変数を用意する。
 これらの変数は、それぞれ 1, 0 のいずれかの値を取るものとする。
 1 は、そのメンバーがそのグループに属していることを意味する。

* [条件1] 各メンバーは、どこかのグループに属する。
 変数の、メンバー行の合計 = 1。

* [条件2] 各グループに属する人数は決まっている。
 変数の、グループ列の合計 = (所定の人数)。

* [条件3] 仲良し同士は同じグループ。
 仲良し同士について、変数をグループ列方向に掛け算した合計 = 1(どこかのグループで1×1となる)。

* [条件4] 嫌い同士は異なるグループ。
 嫌い同士について、変数をグループ列方向に掛け算した合計 = 0(すべてのグループで1×1とはならない)。

* Scalaソースコードはこちら・・・
>> http://brownian.motion.ne.jp/memo/Copris/Groups.scala

好き、嫌いの条件をたくさん入れすぎると、解無しになってしまう。
そうなったとき、1つずつ条件を減らしてゆくと、どこかで解が出てくる。
つまり、誰が我慢すれば丸く収まるのか試すことができる。

■ 環境設定

制約プログラミングのソフトウェアはいくつかあるが、Coprisの良いところは敷居が低いことだと思う。
とにかくScalaさえ動かせれば、Copris自体で覚えるべきことはかなり少ない。
なので、「制約プログラミングとは何ぞや」を知るにはベストなのではないかと思う。
ただ、この「Scalaさえ動かせれば」のところでつまずく人も多いと思うので、Windows上での簡単な導入方法を以下にメモっておく。
(実際「Scala インストール」で検索すると、やれ開発ツールを入れろ、sbtを入れろ、
といった方法がヒットするので、目的に到達する前に息切れしてしまう。スタートはもっと簡単でよい。)

(1). Java runtime version 1.8 以降をインストールする >> http://www.java.com
(もちろんJDKでもかまわない。試しにコマンドプロンプト上で、
  > java -version
 と入力してみて、version 8 以上と出てきたらインストールは不要。)

(2). Scala version 2.11 をインストールする
* Scala Download >> https://www.scala-lang.org/download/
 『CoprisはScala version 2.11で動作する (他のバージョンでは動作しない).』
 とあるので、現在の最新版ではなく、前のバージョンを入手する。
 Scalaダウンロードページの下の方に「Other Releases」「Scala 2.11.12」とあるので、そちらをクリック。
 SCALA 2.11.12ダウンロードページには「DOWNLOAD INTELLIJ」「DOWNLOAD SBT」とあるのだが、
 それぞれ本格的開発向けなので、そこはパスする。
 ページ下の方にある「Other resources」「scala-2.11.12.zip」から、直接zipファイルをダウンロード入手する。

(3). JavaScalaにパスを通す。coprisにクラスパスを通す。
作業するコマンドプロンプト上で、以下のように入力するか、以下のようなバッチファイルを作っておく。
("C:\MyWork"といった箇所は、各人の環境に応じて適切な場所をセットする。)

SET JAVA_HOME=C:\MyWork\java\10
SET SCALA_HOME=C:\MyWork\Scala\scala-2.11.12
PATH=%PATH%;%JAVA_HOME%\bin;%SCALA_HOME%\bin
SET CLASSPATH=.
SET CLASSPATH=%CLASSPATH%;C:\MyWork\Scala\copris-v2-2-8\build\copris-all-v2-2-8.jar

この状況で > scala と入力すると、scala対話プロンプト(REPL)が立ち上がる。

C:\MyWork\Scala>scala
Welcome to Scala 2.11.12 (Java HotSpot(TM) 64-Bit Server VM, Java 10).
Type in expressions for evaluation. Or try :help.
 
scala>


(4). 以降の操作は、Copris入門ページにある通り。
* Coprisによる制約プログラミング入門
>> http://bach.istc.kobe-u.ac.jp/copris/docs/intro-ja.html

チェビシェフの不等式のかんたん理解

どのような標本・確率分布でも・・・平均から 2標準偏差以上離れた値は全体の 1/4 を超えることはなく、
一般にn標準偏差以上離れた値は全体の \frac{1}{n^2} を超えることはない。
    >> wikipedia:チェビシェフの不等式 より.

式で表すと、

  P( |x - \mu| \geq a \sigma) \leq \frac{1}{a^2}
  P() は、カッコの中が成り立つ確率、という意味。
  μは平均。|x-μ| は、個々のデータの値と平均との偏差のこと。
  σ は標準偏差
  a には任意の数を当てはめることができる。
* そんなの常識、あたりまえでない大数の法則 >> http://miku.motion.ne.jp/stories/08_LargeNum.html

このように書くと何だかとても難しいことのように思えますが、実はアタリマエのことを言っているに過ぎません。

● 最も単純な標準偏差1の分布

最も単純な標準偏差1の分布は、データが +1と -1の、2個だけというものでしょう。

 標準偏差σ = √{ (1^2+ (-1)^2) / 2 } = 1。
この状況をチェビシェフの不等式にあてはめると、
『平均0から、1標準偏差以上離れた値は全体の 1/1 を越えることは無い』
つまり、全部のデータを1よりも遠くに引き離すことはできない、ということを言っています。
試しにデータを少しだけ動かして +1.1 と -1.1 にしたならば、それに合わせて標準偏差も 1.1 と大きくなります。
ならば、+1.1 と -0.9 といった具合に動かしてみると、今度は平均が 0.1に上がるだけで、
やはりどちらのデータも標準偏差の1.1を上回る(あるいは-0.9を下回る)ことはありません。


つまり標準偏差とは、データを2個の点で代表させたとき、その広がり方のこと』だったのです。
平均値を『データを1個の点で代表させたとき、その値のこと』だと考えれば、
標準偏差とは、いわば“平均値の2個版”だと見なせます。
データが2個だったとき、チェビシェフの不等式が主張する通り「どのデータも標準偏差を超えることはない」、
・・・そもそも2個のデータの隔たりのことを標準偏差と呼んでいたのだ、と理解できます。

標準偏差が2を越える分布

次に、一部のデータが標準偏差2を越えるような、なるべく単純な分布を考えてみましょう。
2個のデータを +2と -2 に置いて、これらがちょうど標準偏差2に位置するように調整すると、こうなります。

データを +2 と -2 に1個ずつ、あとは0を6個配置する。
最も隔たりの大きい +2, -2 のデータをちょうど標準偏差2の位置に持ってくるには、
標本全体としての標準偏差を1に調整しなければなりません。
それには、±2の広がりを打ち消すだけのデータを平均の0に置く必要があります。
(必ずしも0に置かなくても良いのですが、0に置くのが標準偏差を縮めるには最も効率的です。)
標準偏差を1に保つには、
 { (+2)^2 + (-2)^2 } / (全データの個数) = 1
となるので、(全データの個数) = 8 だと分かります。
このとき、標準偏差2を越える(2以上の)データは8個中2個なので、
確かにチェビシェフの不等式が主張する通り 1/2^2 = 1/4 となっています。

標準偏差がNを越える分布

同じことを、標準偏差3を越える場合で考えると、こうなります。

データを +3 と -3 に1個ずつ、あとは0に16個配置する。
 (全データの個数) = 3^2 × 2
  ・なぜ2乗するかというと、そもそも分散とは各データの偏差の2乗の合計だったからです。
  ・なぜ2倍するかというと、プラス側とマイナス側で2倍になるからです。

標準偏差4を越える場合は、こうなります。

データを +4 と -4 に1個ずつ、あとは0に30個配置する。
 (全データの個数) = 4^2 × 2

標準偏差Nを越えるデータを1個置きたかったなら、N^2 個より多くのデータを0に置く必要がある』
これが、チェビシェフの不等式の意味するところだったのです。

東京→青森、国道4号を通らない山岳ロングライド

東京->青森 732.43km、国道4号を一切通らない山岳ルートを自転車走破!
2018年 9月 23日 AM3:00:00 〜 翌 9月 24日 15:13:43.
時間: 36時間13分43秒 (途中 2時間程度の仮眠)
平均時速: 20.2km
コース: 東京日本橋->江戸川CR->日光->鬼怒川->会津若松->長井->寒河江->新庄->横手->角館->大館->弘前->青森.

ルートラボGPS記録
 (前半)日本橋〜山形新庄 >> https://yahoo.jp/pji0Nv
 (後半)山形新庄〜青森 >> https://yahoo.jp/JX4z0m

■ ことの起こり
今年のシルバーウィークも、どこか遠くに行こうと計画を練っていた。
当初、まだ果たしていない大阪→東京キャノンボールに挑もうと考えていたのだが、
周囲の猛反対に合ってあえなく取り下げとなった。
そこで計画を改めて、できるだけ交通量が少ないルートで遠くまで行くことにした。
(実は同じことを3年前にもやった気がする。。。今後、大阪→東京を走ることはもう無いだろう。)
まっさきに候補に挙がったのは、今年のGWに行った秋田までの再チャレンジだ。
 ・2018/4/28 東京→秋田ロングライド >> d:id:rikunora:20180604
あの時は予想を超える寒さに苦しめられたが、今度は十分な防寒装備で挑んでみようか。
改めて Google Mapを眺めると、秋田までの道をそのまま延長すれば、青森まで達することに気が付いた。
距離を測ってみると、4号線経由と大差無い。
東北の真ん中を突っ切って行くのだから、むしろ距離的には少し短い。
もちろんアップダウンはあるのだが、それでも1000mを越える峠は無い。行けそうだ!

東京→青森は、3年前に4号線経由で行ったことがある。
 ・2015/9/19〜20 東京→青森ロングライド >> d:id:rikunora:20151001
あの時は、とにかく飽きるほど長かったという記憶しか無いのだが、
それでも運良く34時間以内に走りきることができた。
今回は1日半=36時間を目標に青森を目指すことにした。
今回も、着替えと輪行袋をあらかじめ青森中央郵便局に局留めで送っておいた。
ただ、輪行袋については別に緊急脱出用の軽量のものを携行した。
今回はスケジュールぎりぎりの、走行後の翌日午後には出社予定だったので、万が一でも帰れるように備えたわけだ。

■ 機材
・Time VX Elite
・パーツ: Dura-Ace 7800, ワイドレシオFront:50x40, Rear:11x28
・ホイール: Ksyrium SLS
 もうこれ以上はないだろうと思っている、ロングライドの鉄板機材。

・タイヤ: Continental Grand Prix TT
・チューブ: Panaracer R'Air 軽量チューブ
この Grand Prix TT、名前の通り Time Trial用タイヤで、軽くて速い。
全体的に薄いのだが、それでも今までパンク無しの優れものだ。
さらに今回はパンクのリスクを恐れず軽量チューブにしてみた。
この組み合わせの走り心地はすばらしい。

・フロントライト: CatEye Volt300×2 -- バッテリーはVolt400のもの.
リアライト: CatEye Rapid-mini×2
・バーエンドライト: CatEye LOOP2 -- トンネルなどで手を離さずに点灯できる.
Volt300 に Volt400 のバッテリーを付けると、点灯時間が長持ちする。
ライトは前後とも、同じものを2個ずつ装備した。
たとえ1個がトラブったとしても、もう片方で乗り切るためだ。
あと、サイコンは持って行かなかった。(記録はGPS-Watchでとった。)
途中で電池切れになることが分かっていたし、これほどの長距離で、
途中のスピード経過をいちいち気にするのが良くないように思えたので。

■ 服装
ついこのあいだまで猛暑だったのにもかかわらず、東北の夜はかなり冷え込むらしい。
そのため服装についてはかなり悩んだ。
天気予報を参考に、昼は28度〜夜は10度まで耐えられる装備を選んだ。


反射ベスト(Biemme)、防寒用の上着(cannondale, 袖が着脱できてベストにもなる)、ろんぐらいだぁす!ジャージ。


蛍光ウィンドブレーカー(CRAFT)、アンダーシャツ(CRAFT)、タイツ(PearlIzumi)、靴下(PearlIzumi)。


手袋は2重、手が痛くなるので。寒さに備えてのメットインナーキャップは功を奏した。

今回よかった一品は、CRAFTの長袖アンダーシャツ。
下着にしてはかなり高いと思えるお値段だったのだが、実に快適。
下着には贅沢する価値があることを覚えた。
日中はこの長袖アンダーシャツの上に、半袖サイクリングジャージ、反射ベストというスタイル。
夜はウィンドブレーカー、あるいは防寒上着を着用した。

■ 緊急Goods

緊急脱出用の軽量輪行袋(SL-100)。野宿のためのレスキューシート。今回は共に出番が無かった。

■ 経過
・深夜3:00、日本橋道路元標を出発。
 出発時の持参食料、おにぎり×3個、ウィダーインゼリー×2個。
・国道4号を通らないよう、まずは東へ。
 隅田川沿い、スカイツリーの傍らを抜け、国道6号を一路江戸川へ。
 金町から江戸川サイクリングコースに入る。
・江戸川CR、深夜はほばがら空きだが、ごく希に散歩、ランニングに出会う。
 気がつくとハンドルにたくさんの蜘蛛の巣がくっついている。
・境から県道17号に入る。一瞬、雨がぱらついたが、降られることはなかった。
 古河付近の路面は濡れていた。際どいところで雨を回避したようだ。
・県道から細い抜け道ルートで小山へ抜けた。ここでおにぎりを食べる。
・小山郊外「扶桑第一公園」でトイレ休憩、おにぎりを全て食べる。今回、朝からやたらと腹が減る。
壬生町鹿沼を抜け、日光までは微妙な登り。
・日光、大谷川を渡る橋のたもとでトイレ休憩、ウィダーインゼリーを消費。
・鬼怒川バイパス、日塩有料道路を通る。それぞれ50円、20円。
川治温泉トンネル手前で休憩、最後のウィダーインゼリー消費。食料が全て無くなった。
・すばらしい景色の山道、だが何も無い。腹が減ってやばい、ハンガーノックが近い。
中三依(なかみより)温泉「山の幸直売センター」にたどり着き、
 店のおばちゃんに「腹が減って動けない、何か食べ物をください」と訴える。
 親切なおばちゃんが、奥からカボチャの煮物と栗ご飯を出してくれた。
 空きっ腹に染み渡るようにうまかった。おばちゃん、ありがとう!
 しかも無料。申し訳ないので、ジュースにアイス、ビスケット類を購入した。


 (写真は Google Mapより)

・薄曇りから晴に代わり、気温上昇。いいかんじに山王トンネルを越える、標高846m(GPS読み)。
・山間の田んぼの中を下り調子で会津若松へ。途中、下郷付近で自販機休憩。
会津若松では、ちょうど「会津祭り」で賑わっていた。侍のかっこうをした人が練り歩いていた。
 駅前は混雑で通れなかったため、裏手に回る。
・13:20「ファミリーマート会津金川町店」で休憩補給。
・国道121で喜多方市街をバイパスし、「道の駅 喜多の郷」でトイレ休憩。
大峠登り途中、日中ダムで12時間経過。この大峠の道の眺めはけっこう好き、トンネル多いけど。
 だいぶ疲労が溜まってきたが、まだ全行程の1/3と思うと笑いがこみ上げてくる。
・トンネルを越え、山形県道4号、次いで県道8号に入り、川西町へ。
 ここはGWに訪れたとき工事中で走れなかったところ。
 田んぼの中、良い風景の道。このあたりで手持ちのおにぎりを次々と消費。
米坂線今泉駅付近の稲荷神社で休憩。
 自販機ジュースと最後のおにぎり消費。ここではGWのときにも休憩した。
 何気ない場所なのだけれど、ちょうどこの辺りで疲労が溜まるポイントだし、水道もあるし、ここで休みたくなるのだな。
・長井通過、最上川沿いに北上。対面で何台かのロードに出会う。
・荒砥から国道287を最上川沿いの道に入る。ここはGWのときに道を間違えたところだが、今回は大丈夫。
・「道の駅 白鷹ヤナ公園」でトイレ休憩。
 あゆ祭りをやっていたようだが、もう終わりの時間だった。
 ここからライト点灯、ナイトランに入る。
・朝日町、寒河江、河北、国道347に入って村山市、問題なく通過。
 明かりも交通量少ない。この辺りでようやく全行程の半分。
 GWには寒くてたまらなかった記憶があるのだが、今回は気温20〜18度で涼しいくらい。
 走行感覚がまるで違う。
大石田駅前から県道を通って国道13号に入る。
・舟形トンネルあたりで、眠気と疲労が襲ってきた。
 少し先の南新庄駅(無人駅)ホームの待合室で仮眠をとる。
 完全個室、トイレ付き(?!)の快適なところだ。21:00〜30分程度の休憩。


 (写真は Google Mapより)

・22:22「ファミリーマート新庄昭和店」で遅めの夕食。
 おにぎりが喉を通らなくなってきているので、甘いパン食とする。
 (ここで買った予備のおにぎり1個は、結局最後まで食べずに持ち越した。)
・深夜23:30、旧道の主寝坂峠を越える。真っ暗、すれ違うものは何も無い。
・峠を越えてから下り坂、走りやすくなる。
 しかしここから明け方にかけて、疲労と睡眠との戦いとなる。
・0:26「道の駅 おがち 小町の郷」で休憩をとる。
 いよいよ気温が下がってきたので防寒上着を着用。やはりこれが役に立った。
 気温は最も低いところで12度の表示を見た。
 これは1回だけで、あとは16〜18度程度だった。
 (この街道沿いには気温の表示があるのだ。)
 東京の暑さからは想像しにくいところだが、これもGWの経験のおかげ。
・3カ所続く道の駅「おがち 小町の郷」「十文字」「雁の里せんなん 雁太郎」は、
 休憩室が24時間開いていると事前に調査済み。(今回せんなんは通過した)
・次の道の駅「十文字」で眠気が極に達した。
 休憩所に寝られるような場所は無かったのだが、構わずカーペットの床でごろ寝する。
 (GWにやむなく寝た、公園のトイレの中より数倍マシだ。)
 1:30〜2:39、約1:00の仮眠。
・横手の手前で深夜3:00、24時間が経過した。仮眠してもやはり眠い。
 (この感触からすると、24時間東京->秋田は極めて困難だ。)
・美郷から県道11号に入り、角館を目指す。
 この間、眠くて疲れてひたすら長かったという思い出しか無い。
 昼間に走れば、さぞ良いところだったろうに・・・
・角館で夜明けを迎える。
 5:00〜5:34、角館市街を抜けた付近のバス停「元町」で約30分の仮眠をとる。
 雪国のバス停は立派な個室になっているのだ。疲労のあまり熟睡する。


 (写真は Google Mapより)

・目が覚めると辺りは薄曇り。
 ここから北秋田市まで約100km、途中、想像以上に何も無いことに驚く。
 地図を見ると鉄道に沿っているのだが、なぜこんなところを走っているのだろうか。
 疑問に思っていたところ、1両だけの汽車がコトトン、コトトン、と追い抜いていった。
 なんとも感動を覚えた。(秋田内陸縦貫鉄道、知る人ぞ知るローカル線だ。)
・このあたりが最も精神的につらかった。
 上り坂なのでペースが上がらない。
 残り200km(東京-静岡以上!)あると思うと、本当に青森にたどり着けるのかと疑いの念が頭をもたげる。
 幸い輪行袋も持っていることだし、あの電車に乗ってしまえば良いのではないか・・・
・足を引きずるように前進を続けると、いつしか登坂のためチェーン脱着所の前を通過した。峠が近いようだ。
 7:51、やっとのことで大覚野峠を越え、青森まで行く自信を取り戻す。天気も晴れてきた。
・ひたすら山中の長い道を下り、9:40「ローソン 北秋田米内沢諏訪岱店」で補給と休憩。
 ようやく町に戻ってきた。

・鷹の巣から国道7号線に入る。青森まであと105kmの看板を見る。もう少しだ。
 あと少しという気持ちも手伝って、不思議と新たな力が湧いてきた。
 国道7号は交通量が多い。ここまでの山中に比べれば圧倒的。
 (というより、ここまでの交通量が異様に少なかった。)
・大館を通過した先で、大荷物のキャンピング自転車に会った。なんでも神奈川から来たのだとか。
・12:09、矢立峠。これまで越してきた峠に比べれば道も良く、問題なく越せた。いよいよ青森県だ。
・12:28「道の駅 碇ヶ関」で最後の大きな休憩。手持ちの黒糖まんじゅうを食べる。
 (ここは以前、東京->青森ランでも立ち寄ったところだ。立ち寄った時刻もかなり近い。)
・ここから青森までは2通りの道がある。
 (1)国道7号をひたすら走る。(2)県道を通り黒石、青森空港を経る。
 (2)は以前走ったルートなのだが、路面があまり良くなく、青森空港への登りもあることから、今回は(1)を選んだ。
・ところが国道7号は、弘前付近からいよいよ交通量が増え、大型トラックもたくさん通るようになった。
 1カ所、陸橋の階段を上らなければならないところもあった。
 幸い自転車が走れる路肩は広く取ってあったのだが、この路肩はあまり整備されておらず草ぼうぼうで苦労した。
・しかも国道7号にも、青森の手前には山越えがあった。総合的には(2)青森空港ルートの方がよかった。
 これも両方走ってみなければ分からないことなのだが。
・9月24日 15:13:43、ようやく青森県庁前に到着。長い長い道のりであった。
 感想、『もうこれ以上走りたくない(笑)』

■ おまけ
Google がオートバイと認識!


練習のべき乗則

『練習は裏切らない。』
この言葉の真偽に一石を投じる、驚くべき研究結果があります。

* 第50回 練習の効果 >> http://www.pitecan.com/articles/WiredVision/wv50/index.html
木村氏は、吉澤章氏の「創作折り紙」という本で紹介されている「みそさざい」という作品を15万回折り続け、折るのにかかった時間がどのように変化したかを記録しました。

その結果、折るのに要した時間は、練習回数の対数に比例するという、
「練習のべき乗則」(Power law of Practice) に従うことが明らかになりました。
『意外性に魅せられて約10年続けました』なのだそうです。
* 中京大学 人工知能高等研究所ニュースNo.17
>> http://www.iasai.sist.chukyo-u.ac.jp/pdf/iasai_news17.pdf

15万回、10年という努力は並大抵ではありませんが、数十回程度の繰り返しであれば、日常の中でもよく行うことがあります。
そうした作業を繰り返すと、慣れることによって、どれほど作業時間が短くなるのか。試しに測ってみました。
実際に私が試したのは「答案の採点」という作業です。
75枚の答案の採点にかかった時間をストップウォッチで測り、まとめた結果がこのグラフです。

グラフ中に引いた曲線は、(採点の秒数) y = 0.1681 x ^ (-0.269) という累乗曲線。
なるほど、かかった時間はおおむね「練習のべき乗則」に一致しています。
同じ内容ですが、グラフを両対数で描き直したものが下になります。

ここで、直線が「練習のべき乗則」なので、全体の傾向としては合っています。

Wikipediaの「学習曲線」によると >> wikipedia:学習曲線

ピロリとアンダーソンは a 、b の実測値をそれぞれ1.40、0.24と求めた。
  {\displaystyle RT=1.40N^{-0.24}={\frac {1.40}{N^{0.24}}}}
この式はかなり普遍的に成り立つ。

この 0.24 という冪乗の定数 b は、私の測定値では 0.269、なので、確かに近い値になりました。
一方、折り紙の方はグラフから見ると 0.172(あるいは 0.2166)となっており、かなり違っているように見えます。
(係数 a は単位の取り方によるように思えるのだが、、、よくわからん。)
また、以下の論文では b = 0.269 、奇しくも私の値と一致していました。
* スキル学習におけるスランプ発生に対する事例分析的アプローチ
>> https://www.jstage.jst.go.jp/article/tjsai/23/3/23_3_86/_pdf
こうして見ると、べき乗則自体は普遍的に成り立ちそうですが、
その定数値まで普遍的というのは疑念の余地があり、やはり作業内容によって変わってくると思うのです。

■ 気付いたこと

* べき乗則は確かに有用.
たった75回程度でも「練習のべき乗則」が見えてくるとは驚きです。
最初にかなり時間がかかったのは、解答が本当に妥当かどうか、文献にあたって調べた時間などが含まれています。
それが後半になると、似たような解答は覚える、途中の過程まで覚えて一目で点数が分かる、などの効果が実感できました。
練習のべき乗則は、ちょっとした繰り返し作業にも十分有効なのです。

* 時間がかかるのは記述の読み取り.
テストの採点で最も悩むのは、記述式の解答を読み取って、意図を解釈することです。
テスト問題のように限定された状況下でさえ、答は予想以上にバラエティに富んでいます。
「なるほど、そう来たか」と唸るような解答も少なくありません。
こうした解答の意図を汲み取り、解答者がどのように考えたのか推測を巡らせるところが採点の奥深さなのです。
さらに、当たらずとも遠からずといった解答に、どれだけ部分点を配点するかが悩みどころです。
全体として不公平にならないように、こっちに点数を付けたなら、あっちにも点数が付かないとおかしいぞ、
といった調整を図ることになります。これが難しい。
採点時間の上下動は主に、こうした解答の解釈・調整に充てられています。
この点が、折り紙のように均質な作業と、採点のように1つ1つが異なる作業との違いで、採点時間が大きくばらつく理由です。
それでも実際に測ってみると、解釈に悩む時間はせいぜい数十秒程度であることも分かりました。

* 満点と0点は採点が早いか.
良くできた模範答案と、その反対に白紙に近い答案は、採点時間があまりかかりません。
評価に悩むことが無いからです。
正直、みんなが100点取ってくれれば、採点する側はとても楽です。
評価に悩むのは上位層でも下位層でもない、最も数の多い中間層です。
このことは体感的には明らかなのですが、それが採点時間に表れているでしょうか。
そこで、テストの点数と採点時間の関係をプロットしたのが、このグラフです。

中央付近が膨らんでいるように見えなくもないのですが、今ひとつはっきりしません。
それでも「90点以上は採点に2分かかっていない」というのは事実です。
さらに採点への慣れの影響を除くため、(テストの点数)×(べき乗則からの残差)をプロットしてみました。

グラフからはっきりした傾向は読み取れません。
中間層では「大きく時間がかかる場合もある一方、さっさと済んでしまうものもあり、振れ幅が大きい」
というのが事実のようです。

いずれにせよ、もしべき乗則が普遍的なら、最初の伸び方を見て、その後、どの程度練習すればどこまで伸びるかの予測が立つはずです。
この予測をもとに、どこまで練習すべきか、あるいはどこで練習を打ち切るべきかの判断が付くわけで、これは極めて有用なルールと言えるでしょう。

「みそさざい」の折り紙。折ってみました。

指数法則を満たす非連続関数

『さらっと言うと、要は、選択公理を認めると f(mn)=f(m)+f(n) を満たす不連続関数が作れてしまう。』
以下にあった、気になる数学ネタ。
* PRMLガール 〜 文芸部のマネージャーが「パターン認識機械学習」を読んだら
>> http://d.hatena.ne.jp/n_shuyo/20130117/prml
このブログ記事には書籍版があって、「あとがきがわりのACガール」にもう少し詳しい解説があります。

当初は「何のこっちゃ?」と思っていたのですが、最近ようやく意味が分かってきたので、以下につらつらと書いてみます。
PRMLガールでは対数について書かれていましたが、ここでは指数について書きます。)

関数 f が、全ての実数 x, y について、

f(m + n) = f(m)・f(n)

という関係を満たすとき「f は指数法則を満たす」ということにしましょう。
指数法則を満たす連続関数は、いわゆる指数関数 f(x) = a^x しかありません。
なぜかというと、

【Step.1】 自然数 n については・・・
f(n) = f(n-1)・f(1) = f(n-2)・f(1)・f(1) = ・・・ = {f(1)}^n
  ここで a := f(1) と置けば、f(n) = a^n となる。


【Step.2】 0については・・・
f(0) = f(0 + 0) = f(0)・f(0) = f(0)・f(0)・f(0)
  よって、f(0) = 1


【Step.3】 マイナスの数については・・・
1 = f(0) = f(n - n) = f(n)・f(-n)
  よって、f(-n) = 1 / f(n)


【Step.5】 分数については・・・
f(1) = f( 1/n + 1/n + 1/n ・・・ ) = {f(1/n)}^n
  よって、f(1/n) = n√f(1) = n√a
    ※ n が偶数の場合、f(1/n) = ± n√a という正負の2つが考えられるが、
    ※ f が自然に連続となるようなプラスの方を採用することにする。


【Step.6】
  あとは有理数の極限をとって、連続となるように実数全体に拡大する。

なので、滑らかな連続関数に限れば f(x) = a^x となるのですが、
ここで連続関数という制限を取っ払って、病的な非連続関数もありにすれば、
指数法則を満たす関数はもっと他にもあるのではないでしょうか。

そう思って上を見直すと、最後の【Step6】を見直せば別の関数が構成できそうです。
病的な関数を探す(でっち上げる)作戦として、実数全体を、
 ・集合A:有理数と、
 ・集合B:集合Aに含まれない無理数
の2つに分けることを考えます。
そして、
 ・集合Aに含まれる数 x については f(x) = a exp(x)
 ・集合Bに含まれる数 y については f(y) = b exp(y)
のように、A,Bそれぞれに異なる大きさの指数関数をあてがえば、思惑通りの関数が作れそうです。

しかしながら、話はそう単純ではありません。
集合Bに属する無理数の中で、(無理数) + (無理数) = (有理数) となるような数があるからです。
例えば (1-√2)という数も、(1+√2)という数も無理数ですが、
2つを足した (1-√2) + (1+√2) = 2 は有理数となります。
そうなると、集合Bに含まれる数の足し算の結果が集合Aにはみ出してしまい、AとBがきっちり分かれません。

何とかして、集合Bに含まれる無理数同士の足し算を、集合Bの中に閉じ込めておくことはできないか。
そこで登場するのが「同値類」というアイデアです。
実数全体を、足し算の結果がお互いにはみ出さないような、たくさんの(無数の)集合に分けることを考えてみましょう。

実数全体を「有理数×無理数」の巨大な(無限の)表に並べることを考えます。
まず、実数Rに含まれる全ての無理数を横一列に並べます。(数直線を引いて、有理数の点だけを削除します。)
とある1つの無理数 x を選び、もし x の有理数倍 a x が他の無理数 y と一致したら、その一致した先の無理数 y を削除します。
たとえば無理数πを選び出したとき、その有理数倍 2π, 3π, 1/5π, -6/7π・・・ などなどを、ことごとく消せ、ということです。
πという子が1人いれば、他は全て「いらない子」です。

この作業をとことん繰り返せば、ついにはどの無理数有理数倍しても、他の無理数と重ならない状況ができるはずです。
本当に「いる子」だけが残った状況です。
削除の結果、全ての無理数が消えてしまうことはありません。
たとえば π と √2 の有理数倍が重なることは無いので、少なくともこの2つは残るでしょう。
(残った無理数)×(有理数) の巨大な掛け算の表を作れば、その作り方からして、
掛け算の表には(0を除く)全ての実数が、落ちや重なり無く掲載されているはずです。
「掛け算九九」ならぬ、「掛け算有理・無理」というわけです。
※ 同じ考え方で、代数的無理数超越数を掛け合わせた「掛け算代数無理・超越」の表もできると思います。

ここまでサラリと書きましたが、この巨大な掛け算表が本当に作れるかどうかは自明ではありません。
というのも、どの無理数を「いる子」と認め、どれを「いらない子」とするか、決め手が無いからです。
決め手が無い、というのは、いったいどこから手を付けてよいのか、順序立てて処理する手順(アルゴリズム)が無い、ということです。
1つ1つの実数を1列に順序立てて整列することができない、と言い換えても良いでしょう。
たとえば、最大の無理数を「いる子」にしようとしても、最大の無理数というものはありません。
最小の無理数も存在しませんし、最も1に近い無理数もありません。
確かに、πという1個の無理数を取り出してしまえば、π×1だけが「いる子」で、残りは「いらない子」にすることはできます。
有理数を順序立てて処理する手順は存在する。)
しかし、実数の中から全ての無理数を順序立てて取り出す方法が存在しないのです。

『どれも空でないような集合を元とする集合(すなわち、集合の集合)があったときに、それぞれの集合から一つずつ元を選び出して新しい集合を作ることができる』
これを「選択公理」と言います。>> wikipedia:選択公理
上の「有理数×無理数の巨大な掛け算表」は、選択公理を前提として、初めてできることだったのです。
※ 参考: バナッハ・タルスキーのパラドックス >> [id:rikunora:20091021]

※ あと、実数を一列に並べることができるとする「整列可能定理」は選択公理と同値な主張です。(8/16追記)
※ >> wikipedia:選択公理には、選択公理と等価な命題として、整列可能定理、ツォルンの補題、テューキーの補題、などが上げられています。

巨大な掛け算表の、縦一列の並びを「同値類」と言います。
たとえば πの同値類は、2π, 3π, (1/5)π, -(6/7)π・・・ などなど、(有理数)×π という形の数の集まりです。
√2 の同値類は、2√2, 3√2, (1/5)√2, -(6/7)√2・・・ などなど、(有理数)×√2 という形の数の集まりです。
有理数それ自身は、(有理数)×1 という形の同値類に集めることにしましょう。
(この1の列だけは、例外的に1という有理数を持ってきます。1のことを単位元と言います。)
それぞれの同値類には“有理数個”の(加算無限の)要素が含まれています。
そしてこの同値類は、無数に(非加算無限)あります。

・実数 x, y の同値関係を、任意の有理数 a によって x 〜 a y と定義する。
・この同値関係によって、0 を除いた実数 R\{0} を同値類に類別できる。
・(R\{0},・) という乗法群を、有理数の乗法 (Q\{0},・) という部分群で割った剰余群 R/Q が構成できる。

類別を終えた後、とある1つの同値類の中での足し算を考えてみましょう。
たとえば πの同値類の中での足し算、2π+3π=5π, 1/5π+(-(6/7)π)= -23/35π などの答は、
全て同じπの同値類に含まれています。
それというのも、有理数の足し算は有理数の中で閉じているため、
(有理数)×π という形の数同士の足し算の答は (有理数有理数)×π になるからです。
  a π + b π = (a + b) π
どの同値類であっても、同値類の中での足し算は、その同値類の中だけで完結し、他の同値類にはみ出すことはありません。

次に、異なる2つの同値類の間での足し算を考えてみます。
同値類の中には、代表要素(いる子)が π や √2 などといった単純に見えるものと、
(2π+3√2) や (5/6π - 7/8√2) などといった複合的に見えるものがあります。
 ※ (2π+3√2) は (有理数)×πでも、(有理数)×√2 でも無いので、
 ※ πの同値類でもなく、√2の同値類でもない、また別の同値類に入っているのです。

そこで、同値類の代表要素について、次の概念を導入しましょう。

・代表要素 z が、他の代表要素 x1, x2, x3・・・ と、適切な有理数 a1, a2, a3・・・ を用いて
   z = a1 x1 + a2 x2 + a3 x3 + ・・・
 となるとき、z は一次従属であると言う。(複合的に見える)

・代表要素 z が、他の代表要素 x1, x2, x3・・・ と、どんな有理数 a1, a2, a3・・・ を用いても
   z = a1 x1 + a2 x2 + a3 x3 + ・・・
 とはならないとき、z は一次独立であると言う。(単純に見える)

この、一次従属、一次独立という概念によって、同値類の代表要素は2種類に分けることができるはずです。
具体的にどんな操作を行って分類するかは分かりません。非加算無限あるので、想像すらつきません。
でも、分けるための基準ははっきりしているのですから、
とにかく全ての代表要素が一次従属と一次独立に分かれたものとしましょう。
そして、全ての同値類の中から一次従属なものを取り除いて、一次独立なものだけを残しましょう。

こうして残った「縦に有理数×横に一次独立な代表要素」の表を用いて、
(0を除く)全ての実数は、表の中にある数の足し算で一意に表すことができます。
別の言い方をすれば、(0を除く)全ての実数は、表に残った一次独立な代表要素に分解できます。

※ 足し算のやり方が一意であることは、次のようにして分かります。
※ いま、とある数 z が、
※  z = a1 x1 + a2 x2 + a3 x3 + ・・・ と、
※  z = a1' x1 + a2' x2 + a3' x3 + ・・・ のように、
※ 2つの異なるやり方で分解できたとしましょう。(一部の an = 0 であるケースも含めて。)
※ 2つの式の差をとって、
※  (a1 - a1') x1 = (a2' - a2) x2 + (a3' - a3) x3 + ・・・
※ となりますが、これは x1 が一次独立であったという仮定に反します。

一次独立な代表要素のことを「基底」と呼んでいます。
全ての実数は、基底の一次結合((有理数)×(基底)の足し算)として表されます。
幾何学的なイメージを思い浮かべるなら、基底とは、多次元空間の各次元の軸のようなものです。
3次元空間であれば、全ての点は
  ax + by + cz
で表されるでしょう。
(x, y, z) はそれぞれ (縦, 横, 高さ) に相当する基底です。
同じように、実数という空間は無限次元
  a1 x1 + a2 x2 + a3 x3 + ・・・
で表されることになります。
 ※ こうした幾何学的イメージを持てば、ここまで(0を除く)と入れてきた注釈がうまく空間の中に収まります。
 ※ すなわち、空間の原点に0を当てはめればよいのです。

さて、だいぶ目的に近づいてきました。
当初の目的は、指数法則 f(m + n) = f(m)・f(n) を満たす関数作りでした。
実数全体を同値類という無数の集合に分けた後、それぞれの同値類について、適当に異なる定数を加えます。

たとえば f という関数を、以下のように定義します。

・与えられた数 m を、基底 m = a1 x1 + a2 x2 + a3 x3 + ・・・ に分解する。
・各基底ごとに異なる適当な定数 c を定める。
  c はわりと何でもよいのだが、たとえば次のようなルールで定める。
  基底 xn を10進数表記で表し、最初に出てきた 0以外の数字 c(xn) を取り出す。
  たとえば基底がπであれば、π=3.1415…なので、c(π) = 3。
  基底が√2だったら、√2=1.414…なので、c(√2) = 1。
・すべての規定について an xn + c(xn) を足し合わせた結果の指数関数を、関数 f として定義する。
  基底ごとに異なる適当な定数 c(xn) を足す、というところがミソ。
  f(m) = exp( Σ[n]{ an xn + c(xn) } )

例:
 m = 1π + 2√2
 n = 3π + 4√2
のとき、
 f(m + n)
   = exp( 1π + 3 + 2√2 + 1 + 3π + 3 + 4√2 + 1 )
   = exp( 1π + 2√2 )・exp( 3π + 4√2 )・exp(3)・exp(1)・exp(3)・exp(1)
 f(m)・f(n)
   = exp( 1π + 3 + 2√2 + 1)・exp(3π + 3 + 4√2 + 1 )
   = exp( 1π + 2√2 )・exp(3)・exp(1)・exp( 3π + 4√2 )・exp(3)・exp(1)

確かに、この関数は指数法則を満たしているぞ。
この関数 f をグラフに描けば、exp(n)倍という整数刻みで分裂した、穴だらけの指数関数として表されることでしょう。
あるいは、

・適当な定数 c を、基底と同じ値に定める。
 たとえばπという基底に付ける定数はπ。√2という基底に付ける定数は√2。
  f( 1π + 2√2 ) = exp( 1π + π + 2√2 + √2 )

といった方法で関数 f を作ったなら、この関数 f のグラフは(y>0の)平面上をびっしりと埋め尽くす、穴だらけの指数関数となるでしょう。

検索したら、こんな記事があった。なるほど。
* ハメル基底とf(x+y)=f(x)+f(y)をみたす関数
>> http://math-note.xyz/set-theory/application-of-hamel-basis/

フィボナッチ数列と微分方程式の間

フィボナッチ数列の行列表示

フィボナッチ数列とは、1, 1 から出発して、
 1+1=2、1+2=3、2+3=5、3+5=8 のように、
2つの数字を次々と足し合わせてできる数列のこと。
 ※ 0, 1 から出発して、0+1=1、1+1=2、としても、その先は同じになる。
フィボナッチ数列を Fn という記号で表せば、

  F1 = F2 = 1,
  Fn+2 = Fn+1 + Fn  (n = 1, 2, 3, ・・・)

驚きなのは、このフィボナッチ数列の一般項(n番目の数字)が

   Fn = \frac{1}{\sqrt{5}} \left\{ (\frac{1+\sqrt{5}}{2})^n + (\frac{1-\sqrt{5}}{2})^n \right\}

という式になることだ。
なぜこうなるかは、ググれば出てくると思うので、ここでは少し趣を変えて「行列」という見方をしてみよう。

平面上の (a, b) という点を、(b, a + b) という点に移す2×2行列を考える。

   \begin{pmatrix} 0 \quad 1 \\ 1 \quad 1 \end{pmatrix} \begin{pmatrix} a \\ b \end{pmatrix} = \begin{pmatrix} b \\ a+b \end{pmatrix}

この変換行列をAと名付けておこう。
行列Aで、点(1, 1) からスタートして次々に変換を繰り返せば、フィボナッチ数列ができあがる。

   \begin{pmatrix} 0 \quad 1 \\ 1 \quad 1 \end{pmatrix} \begin{pmatrix} 1 \\ 1 \end{pmatrix} = \begin{pmatrix} 1 \\ 2 \end{pmatrix}
   \begin{pmatrix} 0 \quad 1 \\ 1 \quad 1 \end{pmatrix} \begin{pmatrix} 1 \\ 2 \end{pmatrix} = \begin{pmatrix} 2 \\ 3 \end{pmatrix}
   \begin{pmatrix} 0 \quad 1 \\ 1 \quad 1 \end{pmatrix} \begin{pmatrix} 2 \\ 3 \end{pmatrix} = \begin{pmatrix} 3 \\ 5 \end{pmatrix}
    ・・・

変換の様子をグラフにプロットしてみよう。

青い点の折れ線が、(1, 1)からスタートしたフィボナッチ数列だ。
オレンジの折れ線は、試しに(1, 4)からスタートした数列を描いたものだ。
見ての通り、どちらの数列もジグザクを繰り返しつつ、とある一本の直線に漸近している。
実はこの直線の傾きが (1+√5)/2 なのである。この傾きを「黄金比」とも言う。

どのようにして、この直線が出てくるのか。
Aという行列は、
 ・点 (1, 0) を (0, 1) に、
 ・点 (0, 1) を (1, 1) に、
移す変換のことである。
グラフの上では、
 ・赤い実線の矢印を、赤い点線の矢印に、
 ・青い実線の矢印を、青い点線の矢印に、
移す変換となる。(グラフ上では長さ1ではなく、5倍の長さの矢印を描いた。)
この変換によって、ピンク色の正方形は斜め45度で裏返された後、水色の平行四辺形に引き伸ばされた形となる。
この、裏返し->引き伸ばし、裏返し->引き伸ばし・・・が繰り返される様子を想像すれば、
平面全体が、引き伸ばされる方向に流れてゆくイメージが浮かぶだろう。
その流れる先が、黄金比の傾きを持つ直線となるのである。

■ 行列の対角化

この流れる先を計算する方法が「行列の対角化」だ。
対角化とは、ざっくばらんに言えば、行列による変換を回転操作と拡縮操作に分解することだ。
行列Aによる変換を、伸び縮みの方向がちょうど縦横(xy座標軸)となるように、うまい具合に回転する。
(正確に言うと、2本の基底ベクトルをそれぞれ回転して、xy座標軸に重なるようにセットする。
 なので、単なる回転ではなく、ひしゃげた回転とでも言うべきものだ。)
この状態で縦横に(xy座標軸の方向に)拡縮した後、再び逆の回転操作を行って元に戻す。

 (変換A) = (回転P)・(拡縮D)・(逆回転P^-1)


このように、とある変換を回転と拡縮に分解できれば、
変換Aの繰り返しは、以下のように拡縮Dの繰り返しに落とし込むことができる。


  A^n = (P・D・P^-1)・(P・D・P^-1)・(P・D・P^-1)・・・(P・D・P^-1)
    = P・D・(P^-1・P)・D・(P^-1・P)・D・(P・D・P^-1)・・・D・P^-1
    = P・D^n・P^-1

見ての通り、一回前の(逆回転P^-1)と、続く(回転P)が互いにキャンセルし合って、
 (変換A)^n = (回転P)・(拡縮D)^n・(逆回転P^1)
となるわけだ。

拡縮Dをn回繰り返した結果は、少し考えれば分かる。
例えば、横に2倍、縦に3倍拡大する変換をn回繰り返した結果は、
横に 2×2×2・・・ = 2^n倍、縦に3×3×3・・・ = 3^n倍だ。
一般に、横の拡大率をλ1、縦の拡大率をλ2とすれば、

   D^n = \begin{pmatrix} \lambda_1 \quad 0 \\ 0 \quad \lambda_2 \end{pmatrix}^n = \begin{pmatrix} \lambda_1^n \quad 0 \\ 0 \quad \lambda_2^n \end{pmatrix}

※ 縦横の拡縮を行列で書けば対角行列となるので「対角化」と言うわけだ。

問題はどうやって、行列Aをうまい具合に回転と拡縮に分解するか、なのだが、
ここで登場するのが固有ベクトルという概念だ。
固有ベクトルとは、行列の変換によって向きを変えないベクトルのことだ。
上のフィボナッチ数列の例で言えば、点列の流れる先の、黄金比の傾きを持つ直線のことだ。
変換によって向きを変えないベクトルは、それがそのまま平面全体の拡縮の方向を示している。
そのような「向きを変えない不動のベクトル=拡縮方向のベクトル」が見出せれば、
あとはその固有ベクトルをxy座標軸に重なるように持って行くことで、目的を達成できる。

固有ベクトルは必ずしも平面グラフに描けるとは限らないが、複素数まで含めて考えれば、必ずどこかにある。
※ また、複数の固有ベクトルがたまたま重なって一致することもある(重解)。
* 複素数固有ベクトル >> d:id:rikunora:20161208

固有ベクトルの求め方は、やはりググれば出てくるだろうから、簡単に。
固有ベクトルとは、「変換行列Aを施した結果が、やはり元のベクトルのλ倍となっている」ベクトルのことである。

   A \begin{pmatrix} x \\ y \end{pmatrix} = \lambda \begin{pmatrix} x \\ y \end{pmatrix}

0でないベクトル(x, y)に変換を施した結果が0になるためには、変換の行列式が0である必要がある。

   \left\| \begin{pmatrix} a - \lambda \quad b \\ c \quad d - \lambda \end{pmatrix} \right\| = 0
   (a - \lambda)(d - \lambda) - b c = 0   ・・・これを「固有方程式」と言う。

2x2行列の場合、この2次方程式を解くことで、まずλが求まる。
このλのことを固有値と言う。(2次方程式なので、λは2つ出てくる。)
固有値とは、つまり固有ベクトルに沿った拡大率・縮小率のことだ。
得られたλを元の変換に代入して、固有ベクトルが見出せる。
(2つのλに対して、固有ベクトルも2本出てくる。2本の固有ベクトルは互いに直交している。)

こうして得られた固有ベクトルを並べて作った行列が、回転操作Pに相当する。
なぜなら回転操作Pとは、xy座標軸(1,0)と(0,1)を固有ベクトルに重ねる(変換する)操作のことだからだ。
そして、対角行列Dとは、固有値を対角線上に並べた行列になる。
なぜなら対角行列Dは、横にλ1、縦にλ2 拡縮する操作のことだからだ。
最後の逆回転操作 P^-1 は、回転操作Pの逆行列となる。
* 固有ベクトルが直交するのは >> d:id:rikunora:20090307
* 固有ベクトルが直交するのは(2) >> d:id:rikunora:20110203

上の手順でフィボナッチ数列の行列を対角化すれば、結果はこうなる。

   \begin{pmatrix} 0 \quad 1 \\ 1 \quad 1 \end{pmatrix}^n = \begin{pmatrix} - \alpha \quad -\beta \\ 1 \quad 1 \end{pmatrix} \begin{pmatrix} \beta \quad 0 \\ 0 \quad \alpha \end{pmatrix}^n \left\{ \frac{1}{\alpha-\beta} \begin{pmatrix} -1 \quad -\beta \\ 1 \quad \alpha \end{pmatrix} \right\}
  ここで  \alpha = \frac{1+\sqrt{5}}{2}, \quad \beta = \frac{1-\sqrt{5}}{2} .

一見複雑に見えるが、これを整理すると、最初に書いたフィボナッチ数列の一般項の式に一致する。

以上はフィボナッチ数列の話だったのだが、同じ方法が他の漸化式一般にも使えることは明らかだろう。
例えば、3つの数字を次々と足し合わせてできるトリボナッチ数列であれば、次の行列を対角化することで一般項が求められる。

   \begin{pmatrix} 0 \quad 1 \quad 0 \\ 0 \quad 0 \quad 1 \\ 1 \quad 1 \quad 1 \end{pmatrix}  \begin{pmatrix} a \\ b \\ c \end{pmatrix} = \begin{pmatrix} b \\ c \\ a+b+c \end{pmatrix}

 ※ トリボナッチ数列の場合、固有方程式が3次方程式となるので、一般項はかなり複雑になる。
 ※ 数式処理ソフトMathematicaを使って計算した結果を、記事の最後に貼り付けておいた。

■ 漸化式から微分方程式

さて、漸化式をわざわざ行列に書き直したのは、次につなげたかったからだ。
 『漸化式のステップを、うんと細かくしたもは、微分方程式となる。』
漸化式は整数ステップで進行するが、これを極めて小さなステップで進めれば、微分方程式となる。
つまり漸化式の整数nを、実数tに置き換えたものが微分方程式なのである。
具体的に、フィボナッチ数列を細かく刻んだ微分方程式は、こうなる。

  Fn+2 = Fn+1 + Fn    ・・・フィボナッチ数列
  X(t)''= X(t)'+ X(t)   ・・・対応する微分方程式

この微分方程式の解は、以下の通り。

   X(t) = C_1 \, e^{ \frac{\sqrt{5}+1}{2} t } + C_2 \, e^{ - \frac{\sqrt{5}-1}{2} t }

この解は明らかに、フィボナッチ数列の一般項の類似物だ。

フィボナッチ数列と同じように、微分方程式を行列に書き直してみよう。
X(t) の微分 X(t)' を V(x) という記号に置き換える。
すると2階の微分方程式は、このような2本の連立微分方程式に書き直せる。

      V(t) = X'(t)
  X(t) + V(t) = V'(t)

この2本の式をさらに、行列でまとめれば、こうなる。

   \begin{pmatrix} 0 \quad 1 \\ 1 \quad 1 \end{pmatrix} \begin{pmatrix} X(t) \\ V(t) \end{pmatrix} = \frac{d}{dt} \begin{pmatrix} X(t) \\ V(t) \end{pmatrix}

結局のところ、漸化式も(線形)微分方程式も、コアとなる部分は「行列の繰り返し操作」だったのだ。

この微分方程式の解法については、やはりググれば見つかると思うので、繰り返さない。
 ※「定数係数 2階線形同次微分方程式」などでググってみよう。
この微分方程式を解くためには、まず以下のような二次方程式を解く。

   \lambda^2 - \lambda - 1 = 0   ・・・特性方程式

この「微分方程式を代数方程式に置き換えたような式」のことを特性方程式と言うのだが、
この特性方程式の正体は、実は先ほど行列の対角化に見た固有方程式(の同等物)なのである。
コアとなる行列が同じなのだから、解き方も同じなのだ。

特性方程式を解くと、平面全体の拡縮率であるλが得られる。
微分方程式の場合、このλまで求まれば、それで一般解を組み立てることができる。
というのも微分方程式の一般解は、2つの独立な解の組み合わせ(線形結合)でカバーできるからだ。
数列のときに考えた回転操作Pも、所詮は線形結合なので、わざわざ回転させる必要が無い。
どの2つの独立な解を選んでもかまわないのであれば、通常は最も単純な「λ1倍拡縮」と「λ2倍拡縮」の線形結合を選ぶだろう。
 ※ 回転操作は積分定数 C1, C2 が吸収しているわけだ。

なぜ微分方程式の一般解は、線形結合でカバーできるのか。
元を正せば、今行っている操作が全て、ベクトルの足し算と掛け算の組み合わせ(1次式)の範囲内だからだ。
なので、平面であれば、とにかく2本の元になるベクトル(特殊解)さえ見つかれば、あとはその組み合わせで作り出せるわけだ。



■ なぜ  e^\lambda なのか

漸化式と微分方程式を比べたとき、漸化式では、繰り返し操作が対角行列  D^n となっているのに対し、
微分方程式では、特性方程式の解λを  e^{\lambda t} のように、e の累乗とするところが著しく異なっている。
この e (自然対数の底) は、どこから湧いて出てきたのだろうか。

実は、 e という数は、整数ステップを無限小ステップに組み直す過程で立ち現れる数なのだ。
e の素性を知るために、単純な漸化式 An+1 = 2 An を考えてみよう。
1, 2, 4, 8, 16 ・・・ のように、自分自身を2倍にする数列は、
「次の項までの差分が自分自身の値に等しい」という性質を持つ。
  An+1 - An = An
たとえば3番目の項「4」と、その次の「8」の差は、8−4=4となっている。
ただ、これは漸化式が整数ステップだからであって、同じことを小数ステップで行うと結果が違ってくる。

2倍に増える数列を、1年間で100%の利子が付く銀行に例えてみよう。
1年で100%の利子が付く銀行に、1000円の預金を1年間預けたら、2000円になって戻ってくる。
これを半年ステップに分割し、最初の半年で増えた預金を全ていったん引き落とし、再び預け直したら、どうなるか。
 ・まず最初の半年で 1000円×(1+(1/2)) = 1500円に増え、
 ・残りの半年で 1500円×(1+(1/2)) = 2250円となる。
同じように、これを3ヶ月ごとのステップに4分割したならば、
 ・1000×(1+(1/4)) = 1250
 ・1250×(1+(1/4)) = 1562.5
 ・1562.5×(1+(1/4)) = 1953.125
 ・1953.125×(1+(1/4)) = 2441.40625
となる。
1年間をn分割した場合、結果は  \left( 1 + \frac{1}{n} \right)^n 倍に増えることになる。

このnをどこまでも増やした極限  \lim_{n \to \infty} \left( 1 + \frac{1}{n} \right)^n = 2.71828\cdots のことを、我々は e と呼んでいたのである。
e とは、「次の項までの差分が自分自身の値に等しい」数列 An+1 - An = An の、n を無限に小さく刻んだ極限のことだ。
* eについて >> http://miku.motion.ne.jp/precomipo/ (要Flash)

An+1 - An = An に対応する「傾きがその場の値に等しい」微分方程式は、X(t)' = X(t) となる。
この微分方程式の解は  X(t) = C e^t となる。
というより、この微分方程式の答になるような数を探し出してきて、それを e と名付けたのである。

ここまでくれば、なぜ漸化式で  \lambda^n だったものが、
微分方程式 e^{\lambda t} に置き換わるのか、はっきりするのではなかろうか。
整数ステップを無限小ステップに拡張するには、 \lambda^n \left( e^\lambda \right) ^t に置き換える必要がある。
 ※ なお、上の具体例を当てはめると分かるのだが、
 ※ 漸化式の 2^n に対応するものが、微分方程式では (e^1)^t となる。
 ※ つまり漸化式の「λ=2」が、微分方程式の「λ=1」に対応する。
 ※ この違いはどこから来るのかと言えば、微分の X(t)' が、漸化式の差分 An+1 - An に対応するからである。

■ 単振動運動の行列表示

さて、漸化式を細かく刻んで微分方程式ができるなら、逆に微分方程式を細かい漸化式として解釈することも可能となる。
ここでは単振動運動の微分方程式

   X(t)'' = - k X(t)
    (k は振動の固さを表す比例定数、物体の質量 m は1とした)

を例に取り上げてみよう。
この微分方程式を、単純に漸化式に直そうとすると、
  An+2 = - k An
 「現在の値をマイナスにしたものが、2つ後にくる?」
となり、これだけではうまくいかない。
そこで、もう1つ、
 「1つ後の項は、現在の値がそのままコピーされる」※
というルールを付け加えると、漸化式はこのように書ける。
  An+1 = An
  An+2 = - k An
    ※このルールは「慣性の法則」、できるだけ直前の状態を維持することの現れなのだと思う。
この漸化式を行列に直すと、このようになる。

   \begin{pmatrix} 0 \quad 1 \\ -k \quad 0 \end{pmatrix} \begin{pmatrix} X_1 \\ X_2 \end{pmatrix} = \begin{pmatrix} X_2 \\ -k X_1 \end{pmatrix}

この行列を対角化すると、

   \begin{pmatrix} 0 \quad 1 \\ -k \quad 0 \end{pmatrix} = \begin{pmatrix} \frac{i}{\sqrt{k}} \quad -\frac{i}{\sqrt{k}} \\ 1 \quad 1 \end{pmatrix} \begin{pmatrix} -\sqrt{k}i \quad 0 \\ 0 \quad \sqrt{k}i \end{pmatrix} \begin{pmatrix} -\frac{\sqrt{k}i}{2} \quad \frac{1}{2} \\ \frac{\sqrt{k}i}{2} \quad \frac{1}{2} \end{pmatrix}

いちいち  \sqrt{k} が出てくるので、最初から  k = \omega^2 と置いておけば、結果がすっきりする。

   X(t)'' = -\omega^2 X(t)

   \begin{pmatrix} 0 \quad 1 \\ -\omega^2 \quad 0 \end{pmatrix} = \begin{pmatrix} \frac{i}{\omega} \quad -\frac{i}{\omega} \\ 1 \quad 1 \end{pmatrix} \begin{pmatrix} -\omega i \quad 0 \\ 0 \quad \omega i \end{pmatrix} \begin{pmatrix} -\frac{\omega i}{2} \quad \frac{1}{2} \\ \frac{\omega i}{2} \quad \frac{1}{2} \end{pmatrix}

この微分方程式の解は、対角行列にある固有値を e の肩に乗せたものの、線形結合となる。

   X(t) = C_1 \, e^{ - \omega i t } + C_2 \, e^{ \omega i t }

固有値虚数の場合、固有ベクトルは実数の平面上には描けない。
その場合、解X(t)の流れは回転となり、結果は正弦波 Sin, Cos の組み合わせとなる。


そして、このように行列に分解してみると、なぜ単振動の角速度が ω^2 のように二乗となるのか、理由が見えてくる。
そもそも2階微分とは「小さく2回変換を繰り返す」ことだったのだから、
1回目の変換でω、2回目の変換でω^2となるわけだ。

力学では、位置Xの2階微分が力に等しい(ニュートンの運動法則)とされている。
これを数列としてとらえるなら、力とは「隣り合う3つの項の差分」ということになる。
  Force = (Xn+2 - Xn+1) - (Xn+1 - Xn)
図に描くと、こういうことだ。

力は、傾きが変化するその点にかかる。

この考え方をもとに、単振動のメカニズムを図解すると、こうなる。


オイラーの公式の行列としての意味

ところで、なぜ固有値虚数の場合は回転となるのだろうか。
理由は「オイラーの公式」でググれば出てくると思うが、ここでも行列流の解釈をしてみよう。

オイラーの公式」とは、 e という数が無限小ステップで立ち現れたのと同じ考え方を、2×2行列にあてはめたものなのである。
90度回転  \begin{pmatrix} 0 \quad -1 \\ 1 \quad 0 \end{pmatrix} を対角化すると、

   \begin{pmatrix} 0 \quad -1 \\ 1 \quad 0 \end{pmatrix} = \begin{pmatrix} -i \quad i \\ 1 \quad 1 \end{pmatrix} \begin{pmatrix} -i \quad 0 \\ 0 \quad i \end{pmatrix} \begin{pmatrix} -\frac{i}{2} \quad \frac{1}{2} \\ \frac{-i}{2} \quad \frac{1}{2} \end{pmatrix}

90度回転を無限小ステップに刻んだものは、これまでと同じ考え方で、対角行列に現れた固有値を e の肩の上に乗せたものになるだろう。

   \begin{pmatrix} 0 \quad -1 \\ 1 \quad 0 \end{pmatrix}^t = \begin{pmatrix} -i \quad i \\ 1 \quad 1 \end{pmatrix} \begin{pmatrix} e^{-i t} \quad 0 \\ 0 \quad e^{i t} \end{pmatrix} \begin{pmatrix} -\frac{i}{2} \quad \frac{1}{2} \\ \frac{-i}{2} \quad \frac{1}{2} \end{pmatrix}

この行列の繰り返し操作をまとめると(そのまま掛け算すると)、こうなる。

   \frac{1}{2} \begin{pmatrix} e^{-i t} + e^{i t} \quad -i ( e^{-i t} - e^{i t} ) \\ i ( e^{-i t} - e^{i t} ) \quad e^{-i t} + e^{i t} \end{pmatrix}

この行列は「90度回転の無限小ステップ」という意味から考えて、小数の角度tの回転変換に一致するはずだろう。

   \begin{pmatrix} cos(t) \quad -sin(t) \\ sin(t) \quad cos(t) \end{pmatrix}

2つの行列の要素を比較して、以下のオイラーの公式を得る。

   cos(t) = \frac{1}{2} ( { e^{-i t} + e^{i t} } )
   sin(t) = \frac{i}{2} ( { e^{-i t} - e^{i t} } )
   e^{i t} = cos(t) + i sin(t)


■ 減衰振動運動の行列表示

単振動が分かったついでに、さらに、抵抗のある減衰振動についても行列化してみよう。

   X(t)'' = -\omega^2 X(t) - 2 r X'(t)
      (ωは角速度、r は抵抗の大きさを表す係数)

この微分方程式を行列に直すと、

   \begin{pmatrix} 0 \quad 1 \\ -\omega^2 \quad -2r \end{pmatrix} \begin{pmatrix} X_1 \\ X_2 \end{pmatrix} = \begin{pmatrix} X_2 \\ -\omega^2 X_1 - 2r X_2 \end{pmatrix}

対角化すると、

   \begin{pmatrix} 0 \quad 1 \\ -\omega^2 \quad -2r \end{pmatrix} = \begin{pmatrix} \frac{ -r + \sqrt{r^2-\omega^2} }{\omega^2} \quad -\frac{ +r + \sqrt{r^2-\omega^2} }{\omega^2} \\ 1 \quad 1 \end{pmatrix} \begin{pmatrix} -r + \sqrt{r^2-\omega^2} \quad 0 \\ 0 \quad -r + \sqrt{r^2-\omega^2} \end{pmatrix} \begin{pmatrix} \frac{\omega^2}{2 \sqrt{r^2-\omega^2}} \quad \frac{1}{2} (1 + \frac{r}{ \sqrt{r^2-\omega^2} }) \\ - \frac{\omega^2}{2 \sqrt{r^2-\omega^2}} \quad  \frac{1}{2} (1 - \frac{r}{ \sqrt{r^2-\omega^2} })  \end{pmatrix}

この微分方程式の解は、対角行列にある固有値を e の肩に乗せたものの、線形結合となる。

   X(t) = C_1 \, e^{ -r - \sqrt{r^2-\omega^2} t } + C_2 \, e^{ -r + \sqrt{r^2-\omega^2} t }

ここでも抵抗の係数 r を、なぜ最初から 2r と置くとすっきりするのか、理由が分かってくる。
抵抗は2回の変換について、1回目で r、2回目でもう +r だけかかる。
なので、合わせて 2r となるわけだ。
(ωが掛け算でかかっていたのに対し、抵抗rは足し算でかかる。)

■ まとめ

・漸化式のステップを、うんと細かくしたもは、微分方程式となる。
・漸化式も(線形)微分方程式も、コアとなる部分は「行列の繰り返し操作」
・整数ステップを無限小ステップに拡張するには、 \lambda^n e^{\lambda t} に置き換える。

数列、行列、微分方程式、この3つは合わさって1つの「線形」という世界を形作っている。
というより、もともと世界は1つなので、3者は1つの世界を別の側面から切り取った姿なのである。
大学や専門学校では1つ1つを学ぶのに息切れして、なかなか3つ合わせた1つの世界までたどり着けない。
それでも線形世界は、1つだけ、あるいはバラバラなまま終わってしまうには余りにも惜しいほどの視点を提供する。
3つは1つなのだと心してほしい。

■ おまけ:トリボナッチ数列の対角化


なぜx2乗の微分は2xなのか

x2乗と聞いて、こんな放物線を思い浮かべたり、

微分と聞いて、こんなイメージが浮かぶ人は、かなり数学を勉強したのだろうと思う。


でも、ここではいったん放物線のことは忘れて、
x2乗とは、以下のように一辺xの正方形を順番に並べたものだと考えてみよう。

ここで、とある正方形と、すぐ隣の正方形との差分は、こんな風になる。

うんと近くの隣同士だったなら、コーナーにある小さな h^2 の正方形は無視できて、
実質的な差分は、上辺と右辺に張り付いている細長い「2つの、長さxの長方形」になる。
差分は 2xh なのだが、これを「すぐ隣までの間隔h」で割ると、2xだけが残る。
だから、
  『x^2 の、すぐ隣同士の差分は 2x』
これが x2乗の微分
  (x^2)' = 2 x
の意味だ。

同じように、x3乗という関数を、こんな風に一辺xの立方体を順番に並べたものだと考えてみよう。

この立方体の列で、すぐ隣の立方体との差分は、こんな風になる。

うんと近くの隣同士だったなら、コーナーにある小さな h^3 の立方体と、
辺に張り付いている3本の細長い棒のような立体 xh^2 は無視できて、
実質的な差分は、上面と右面と背面に張り付いている薄い3枚の正方形、「3つのx^2」になる。
差分は 3hx^2 なのだが、これを「すぐ隣までの間隔h」で割ると、3x^2だけが残る。
だから、
  『x^3 の、すぐ隣同士の差分は 3x^2』
これが x3乗の微分
  (x^3)' = 3 x^2
の意味だ。

同じように、x4乗という関数は、4次元の立方体を順番に並べたものだと考えられる。
4次元の図形は直接描けないけれど、4次元立方体の隣同士の差分を想像すると、
4つの次元それぞれの面(体と言うべきか)に張り付いた
(4次元世界から見れば)薄い4個の3次元立方体となるはずだ。
なので、
  (x^4)' = 4 x^3

この関係は次元がいくつ増えようとも変わらない。
n次元の立方体の差分とは、n次元立方体のn個の表面に、(n−1)次元の薄い立体が張り付いたものだ。
なので、
  (x^n)' = n x^(n-1)
これが、べき乗の微分の公式だ。

・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・
いま、上で見てきたことの逆が「積分」だ。
微分「すぐ隣同士の引き算」なら、積分とは「すぐ隣同士の足し算」のことだ。
x^2 の積分が (1/3) x^3 になるとは、
正方形を重ね合わせて作った立体の体積が、一辺 x の立方体の 1/3 になる、ということだ。
つまりこれは、四角錐の体積のことだ。


1つ次元を下げて、
x の積分が (1/2) x^2 になるとは、
長さ x の棒を、短い方から順に並べて作った図形の面積が、一辺 x の正方形の 1/2 になる、ということだ。
要は三角形の面積のことだ。
(※ただし、ここで x は0から数え始めているので、積分定数Cというものを省いている。)

ところで、四角錐の体積は、なぜ (1/3) になるのだろうか。
イメージしやすいのは、こんな風に6個のピラミッド型を組み合わせた形だろう。

ピラミッドの高さは、立方体全体の半分になっているので、
 (四角錐の体積) = (ピラミッドの体積)×2 = (立方体の体積)÷6×2 = (立方体の体積)×(1/3) = x^3・(1/3)

このピラミッドのイメージは、さらに次元を上げても同じことだ。
4次元の世界で、底面ならぬ“底立体”が3次元立方体となっている“4次元ピラミッド”を、
4つの各次元ごとに2個ずつ8個(4×2=8)組み合わせると、4次元の立方体 x^4 ができあがる。
ということは、4次元の(超)四角錐の体積は、4次元立方体の(1/4)になっている。
以降、n次元世界の(超)四角錐の体積は、n次元立方体の体積の (1/n) となる。
  ∫x^n dx = (1/(n+1)) x^(n+1)

さて、いま“ピラミッド型”の立体をイメージしたのだが、
もう少し想像力を働かせると、四角錐を直接組み合わせて立方体を作る様子をイメージできる。
いま一度、3次元の四角錐に戻ろう。
正方形を並べて作った ∫x^2 dx の四角錐を、3個組み合わせて、1つの立方体を形作る様子がイメージできるだろうか。
白状すると、私自身はこのイメージを頭の中だけで描くことができず、ボール紙で立体を作って確かめた。

立方体を対角線に垂直な面で切ると、切り口は正三角形(場所によっては正6角形)となる。
その正三角形の各辺が、それぞれ1個の四角錐に対応する、といえば分かりやすいだろうか。
よく分からなければ、ボール紙か粘土か何かで、とにかく立体を作ってみれば分かる。
こういうときに常々思うのは、人間(あるいは凡人)は3次元の形をはっきりとは把握していなくて、
分かるのはせいぜい2.5次元程度なのではないか、ということである。

より高次元の積分は、この「四角錐->立方体」イメージの延長上にある。
4次元であれば、3次元の立方体を積み重ねて作った超四角錐を、
4つの方向から4個組み合わせて、4次元立方体を作ることができる。
n次元であれば、(n−1)次元の超立方体を積み重ねて作った超四角錐を、
n個の直交する軸の方向から組み合わせて、n次元立方体を作ることができる。


・・・といった内容のことを、この本に著しました。

実用のための「微積」と「ラグランジアン」 (I・O BOOKS)

実用のための「微積」と「ラグランジアン」 (I・O BOOKS)


* おまけ * Σのイメージ *
「3つの四角錐で立方体を作る」様子がイメージできたなら、
その応用として、Σn^2 = 1^2 + 2^2 + 3^2 + 4^2 + ・・・が想像できるだろうか。
これが、Σn = 1 + 2 + 3 + 4 +・・・であったなら、イメージは難しくない。

要は三角形の“ブロック版”だ。
Σn^2 とは、四角いブロックを積み上げて作った階段状の立体となる。
以下が苦労の末描いた図なのだが、これで分かるだろうか。
(さすがにこれはボール紙で作る訳にはいかなかった・・・
 7/10: その後、積み木で立体を作ることができた。下の追記参照。)

階段状の立体を3個組み合わせると、n×(n+1)×(n+1)の直方体ができるのだが、
その直方体の1枚の表面は、平面の三角形の階段が欠けた形となっている。
この形から、
 Σn^2 = { n(n+1)^2 - ( n(n+1) / 2 ) } / 3
    = n(n+1)(2n+1) / 6
という高校で教わるような公式が得られる。

さらに、この直方体の塊を2つ組み合わせると(つまり階段状の立体を6個組み合わせると)
「三角形の階段が欠けた部分が」互いに噛み合わさって、n(n+1)(2n+1)の細長い直方体ができる。
この式にある(2n+1)の+1は、噛み合わさった面の1列を表している。
なので結局、Σn^2 = n(n+1)(2n+1) / 6 となるわけだ。

そんな階段状の立体の組み合わせでも簡単にイメージできるという想像力豊かな方は、
4次元の世界で、4つの階段状の立体が組み合わさる様子を考えてはいかがだろうか。
私はここで挫折したが。
 Σn^3 = n^2 (n+1)^2
という結果から推測するに、4次元の直方体では表面が欠けることなく綺麗に収まるらしい。
5次元より先は、もはや直観の及ばない世界と言うしかない。

* 積み木で作るΣn^2 (7/10追記) *

100円ショップで売っていた積み木ブロックで、Σn^2 を作ることができた。
まず、n=3 の塊を3個作る。

それらを組み合わせると、こんな形になる。

これを2つ作ってみると、2つの塊が合体できることに気付く。

合体させると、n・(n+1)・(2n+1) の直方体ができる。

これが、Σn^2 = n(n+1)(2n+1) / 6 の意味だ。
(3枚目の写真と色が食い違っています。一度バラバラになってしまったので、再度やり直した。
実際に積んでみると、けっこう難しい!)