「フランクリンの凧問題」というらしい

ここのところ記事を更新していなかったけど、実は、以前書いた三角形の記事の「予想」は正しいと、H.Nakao さんという方からご指摘をいただいて、コメント欄でそれをずっと教えてもらっていました。
ただ結局、コメントにある「Z 係数」のあたりはどうしても分からなかったので、その部分は自分なりに書き換えて、後日ぼくも自分で証明を付けてみようと思っています。
…といっても、近いうちに載せるであろう「証明」は Nakao さんのものであることには変わりありません。
毎日毎日、質問に付き合っていただいて、ありがとうございました。感謝いたします。


主張をもう一度書いておきます↓。元記事と同じです。
※ それと、以後、表記は元記事と同じとします。

主張

下図において、 \triangle ABC AB=AC二等辺三角形とし、0 < x,  y < 80 とする。
このとき、 x, y, z 共に整数の組は、自明なもの *1 を除くと、以下の 16 通りである。

  1. x = 50, y = 20, z = 60
  2. x = 65, y = 25, z = 85
  3. x = 60, y = 30, z = 80
  4. x = 50, y = 40, z = 60
  5. x = 20, y = 50, z = 10
  6. x = 40, y = 50, z = 30
  7. x = 60, y = 50, z = 80
  8. x = 70, y = 50, z = 110
  9. x = 30, y = 60, z = 10
  10. x = 50, y = 60, z = 30
  11. x = 65, y = 60, z = 85
  12. x = 70, y = 60, z = 110
  13. x = 25, y = 65, z = 5
  14. x = 60, y = 65, z = 40
  15. x = 50, y = 70, z = 10
  16. x = 60, y = 70, z = 20


証明のポイント

ぼくが第一に気を付けてるのは、証明に用いるプログラムの中に、「小数点値の等号比較」は絶対に入れないという点です。
元記事は何故×だったかというと、

3. z = 1, 2, ... , 179 において、f(z) = 0 となる z を探す。

の部分を数値計算に落とし込むときに、double 型の等号比較を使ったから。
小数点値の等号比較が入った時点で、それは「近似操作」であり、数学の証明としては不適当なんですよね。


やろうとしている証明は、大まかに言うとこんな感じです。

  • すべての x, y = 1, 2, ... , 79、z = 1, 2, ... , 179 について、以下の操作を行う。
    1. 各 x、y、z に対して定まる f :=「  \frac{d-e}{d-b} \cdot \exp(2 \pi i \cdot \frac{z}{360})虚数部分」*2 を、ζ:=  \exp(\frac{2 \pi i}{360}) の Q 係数 *3 有理式として表す。
    2. 「Q 係数有理式 = 0」を変形して、同値な関係式「ζの Q 係数多項式 = 0」とする*4。この「ζの Q 係数多項式」を g(ζ) (∈Q[ζ]) とする。
    3. ζの Q 上最小多項式 T^96 + T^84 - T^60 - T^48 - T^36 + T^12 + 1 を φ(T) と書く *5 。また、g(T) (∈ Q[T]) を φ(T) で割った余りを ψ(T) (∈ Q[T]) とする。ここで、「 g(ζ) = 0 ⇔ ψ(T) = 0 *6」なので、「g(ζ) = 0」を調べるのには「ψ(T) が 0 である」ことを調べれば良い。

元記事で double 値の等号比較をしてたところを、有理数係数多項式の係数が 0 かどうか、つまり、有限個の有理数列の等号比較に帰着させているわけです。
(有理数は整数 2 つの組を "ある同値関係" で割って定義されるので、小数点値を出さずに済むはずです。)
数学ではよく使ってた手だけど、こうやって数値計算上の困難を解決するのに使うのは思いつきませんでした。


証明に用いるプログラムは

  • Q 係数多項式が用意されている
  • 多項式の剰余」が用意されている

ということから Maxima *7 で作るつもりです。
ただ、

  • Maxima の内部で、本当に小数点値を出さずに計算しているのか?
  • Maxima が使っている定義と、ぼくが使っている定義とはホントに同じなのか?(多項式が「零である」ことの定義とか)

という懸念 (※) は未だ残っているんですけどね。。

感想

まだ完全に解いたわけじゃないけど、とりあえず現時点で思ったことをメモしておきます。

  1. いろんな例を計算するのに使うならともかく、証明の部分にコンピューターを使うのはなかなか厄介だと思った。(※)の部分のように、「この言語 (もしくは数式処理システム) は本当に自分の証明に使えるのか?」という部分のチェックが本当に面倒そう。
  2. 言語 (もしくは数式処理システム) バージョンが上がるごとに、1 のチェックをしなくてはいけないのだろうか??正直それは避けたいなぁ…。
  3. 中学の知識で解けるはず?の問題が、こんな代数的な問題に変形されてしまった。。でも、個人的には、極一部の人にしか見つけられないパズル的な、アクロバティックな解法より、(同じ『極一部』であっても) 「お作法」を身につけさえすれば自然に見えてくる解法のほうが美しいと思う。


*1:x=y=z のもの。

*2:元記事と表記を変えています。

*3:有理数係数」のことです。

*4:上手く説明できなくて申し訳ないんだけど、整数 a, b, c, d (a, c ≠ 0) における「b/a + d/c = 0 ⇔ (bc + ad)/ac = 0 ⇔ bc + ad = 0」と同じ操作です。

*5:この値は Nakao さんの計算結果をそのまま引用させていただいています。

*6:この部分は大学の数学の知識を使っています。

*7:自分の環境では Maxima-5.15.0 です。

「espressione」を最初に音楽用語として使ったのは、モーツァルトなのか??

「espressione」は、「白水社 新伊和辞典」によると、

(女性名詞)

  1. 表現、表現法
  2. 語句、語法
  3. 表情豊かな顔つき
  4. (数学) 式

という意味で *1、英語でいうところの「expression」に対応する言葉だとぼくは思っています。
最近知ったのですが、この言葉は、割と新しい言葉らしく、今道友信さんの「美について」

美について (講談社現代新書)

美について (講談社現代新書)

という本には、こう書かれていました。

このようにして、内的な世界の自己提示として、自己から絞り出すようにしてその感懐を物質の世界に投影するのが、芸術創作の新しい理念として登場してきた "表現" である。この "表現" に対応する言葉は、周知の通りに、古典ギリシア語や古典ラテン語にはなく、エクスプレシオ (expressio) というラテン語である。その当初の意味は果物の果汁を絞り出すという農業用語であり、意外に思う人もいるかと思うが、ゲーテの時代にも、ドイツ語ではまだ "表現" にあたるアウスドルック (Ausdruck) という言葉すらなかった。"表現" は西洋の芸術思想に関する限り、十九世紀から二十世紀にかけてはじめて美学上の市民権を獲得した現代的理念に過ぎない。*2

そうすると、当時、現代的であった「espressione」という単語を最初に音楽用語に使ったのは誰なのか?という疑問が自然に湧いてくるのですが、現時点では、モーツァルトピアノソナタ K.310 の第 2 楽章 (Andante cantabile con espressione)

で使ったのが最初なのではないか…?と推測しています。
それにしても、、何度聴いてもいい演奏ですよね…*3


あと不思議なことに、モーツァルトが「espressione」という単語を使ったのは、現在調べた範囲では K.310 だけのようです (※) *4
K.310 が作曲された 1778 年夏は、モーツァルトの母アンナの亡くなった時期なんだけど、それとは関係があるのでしょうか??
そのちょっと後の世代の作曲家であるベートーヴェンが、「espressione」をしょっちゅう使ってた *5 のとは対照的です。


現時点で調べたのは、K.310 以前ではこの 5 冊

だけなので、もっと調べて信憑性を高めていかねば。
イタリア、フランスの音楽について全然調べてない (うちに楽譜がない…) ので、リサーチ不足と言われても仕方ないですね。。

メモ

(※) の主張の正否については、モーツァルトの全楽譜が公開されているわけだし、
参考:「NMA オンライン」 redhost24-001.com
時間さえかければきちんと調べることが出来るだろうけど、さすがに大変だよなあ。。
楽譜の画像情報から文字情報を検索する、みたいな技術ってないのかなぁ。

*1:「表現」と「式」を同じ言葉で表現する、ってのは、日本人のぼくから見るとどうも不思議な感じがする。。

*2:p85-p86 より引用。ちなみに、「""」は、ぼくが勝手に付けました。

*3:上の動画はリパッティの演奏です。

*4:「cantabile」とか「sotto voce」とか、歌に由来すると思われる単語はよく使ってるんですけどね。

*5:ソナタ第 1 番で、既に使っています。

予想

中学数学に大苦戦! - ひらのの日記
の続きです。


前の記事の問題 1 で不必要に時間をかけすぎたせいか、ただ 1 問解いただけではなんだか報われない気がするので、何か適当な結果を捏ち上げるために (笑) こんな予想を立ててみました。

予想

下図において、 \triangle ABC AB=AC二等辺三角形とし、0 < x,  y < 80 とする。
このとき、 x, y, z 共に整数の組は、自明なもの *1 を除くと、以下の 16 通りである。

  1. x = 50, y = 20, z = 60
  2. x = 65, y = 25, z = 85
  3. x = 60, y = 30, z = 80
  4. x = 50, y = 40, z = 60
  5. x = 20, y = 50, z = 10
  6. x = 40, y = 50, z = 30
  7. x = 60, y = 50, z = 80
  8. x = 70, y = 50, z = 110
  9. x = 30, y = 60, z = 10
  10. x = 50, y = 60, z = 30
  11. x = 65, y = 60, z = 85
  12. x = 70, y = 60, z = 110
  13. x = 25, y = 65, z = 5
  14. x = 60, y = 65, z = 40
  15. x = 50, y = 70, z = 10
  16. x = 60, y = 70, z = 20


x と y をひっくり返しただけ、というのもあるので、この条件下で、本質的に異なる非自明な (x, y, z) の組はたった 8 組だけ、 という予想になります。
もしこれが正しいとすると、、中学で出てきた図形の角度を求める問題って、かなり特殊な状況下での出来事だったんですね。
中学の図形問題を作るのって、大変だったんだなぁ。。

メモ

この「予想」はどうやって立てたかというと、アイデアは前の記事と全く同じで、あとは C++ でコードを書いて、それを実行して出しました。
具体的には、以下の通りです *2


アルゴリズム

  1.  \frac{d-e}{d-b} を x、y で表す。
  2.  \frac{d-e}{d-b} \cdot \exp(2 \pi i \cdot \frac{z}{360})虚数部分は z の関数になるが、それを f(z) と書く。
  3. z = 1, 2, ... , 179 において、f(z) = 0 となる z を探す。
  4. 1 〜 3 の作業を、すべての x, y = 1, 2, ... , 79 について行う。



この 8 通りについては、前の記事と同じように証明すれば良いと思うんだけど (えげつなさそうなのも幾つかあるけど)、問題は、

  • 数値計算で f(z) ≠ 0 と結論付けてしまったものが、ホントに f(z) ≠ 0 なのか?

というところだと思います。
C++ のコードのほうでは、「f(z) の絶対値が DBL_EPSILON より小さいかどうか」でアルゴリズムの 3 の部分の判定を行っているのですが、

  • 数値計算の結果、絶対値が DBL_EPSILON より大きくなってしまったけど
  • 実際は 0 に等しい

ということもありうるわけですからね *3
その辺りのチェックが難しそうです。


もし何か進展があったり、間違いを見つけたりしたら、更新していくことにします。

*1:x=y=z のもの。

*2:表記は前の記事と同じものとします。

*3:ぼくは数値計算についてはド素人なので、適当なことを言ってる可能性が高いです。

中学数学に大苦戦!

問題 1

下図において、 \triangle ABC AB=AC二等辺三角形とする。
このとき、x を求めよ。

この手の問題なら、複素数平面で考えて偏角出せばカンタンに解けるだろう、と思ってたのですが、、
計算してもしても式が簡単にならない!
sin(20^{\circ}) の 10 何次式、とか出てきてましたからねぇ。。

1 週間近くずーーっと、ひらのの日記 の準備も、 5 月に弾くパルティータ第 2 番の譜読みもせずに、この計算ばかりしてたのですが、全然進みませんでした。


一昨日、ついにこの状況に我慢ができなくなって、分度器を買ってきてそれらしい図を描いて、実際に角度を測ってみました。
そしたら 16 度っぽい数字になったので、この辺の値にならないかな〜と再度計算したけど、上手く行かず。。
そんなとき、たまたまこのサイト↓
解答編
を見つけて愕然としました。
どうやら、問題 1 の設定が違ってたようです。がーん…。

  •  BCE = 50^{\circ}

のときに、

  • x = 30

になるみたいです。

この証明 (∠ BCE = 50^{\circ} の場合の証明) を以下に付けておきます。
おヒマな方はどうぞ。
あと、間違いなどがあったら教えてください。(急いで書いたので。。)


その前に、表記と方針を書いておきます。

表記

  1. BC=1 とする *1
  2. ABC\ldots に対応する複素数を、それぞれ、abc\ldots とする。
  3. c-b=1 とする。(BC が実軸と平行、という意味)


方針

  1. \frac{d-e}{d-b}偏角を求める。
  2. その為に、辺 BD の長さ、辺 BE の長さを、初等幾何の方法で求める。


証明:ステップ 1 (辺 BD の長さ、辺 BE の長さを求める)

\triangle BCE は、BC=BE二等辺三角形になるので、
BE=1
になる。
\triangle BCD について正弦定理を適用すると、
 \frac{BC}{\sin(2 \pi \cdot \frac{40}{360})} = \frac{BD}{\sin(2 \pi \cdot \frac{80}{360})}
より、
 BD = 2 \cos(2 \pi \cdot \frac{40}{360})
となる。

証明:ステップ 2 (複素数平面の問題として証明する)

まず、ステップ 1 より、

  • e-b=\exp(2 \pi i \cdot \frac{80}{360})
  • d-b= 2 \cos(2 \pi \cdot \frac{40}{360}) \cdot \exp(2 \pi i \cdot \frac{60}{360})

である。
よって、x = 30 を示すためには、
 \frac{d-e}{d-b} = \frac{(d-b)-(e-b)}{d-b} = 1 - \frac{1}{2 \cos(2 \pi \cdot \frac{40}{360})} \cdot \exp(2 \pi i \cdot \frac{20}{360})
が、
 r \cdot \exp(2 \pi i \cdot \frac{-30}{360}) (r は実数)
という形で書けることを示せばよい *2
つまり、
 (1 - \frac{1}{2 \cos(2 \pi \cdot \frac{40}{360})} \cdot \exp(2 \pi i \cdot \frac{20}{360})) \cdot \frac{1}{\exp(2 \pi i \cdot \frac{-30}{360})} =   \exp(2 \pi i \cdot \frac{30}{360}) - \frac{1}{2 \cos(2 \pi \cdot \frac{40}{360})} \cdot \exp(2 \pi i \cdot \frac{50}{360}) …(★)
が実数、言い換えれば虚部が 0 であることを示せば十分。


(★) の虚部は
 \sin(2 \pi \cdot \frac{30}{360}) - \frac{1}{2 \cos(2 \pi \cdot \frac{40}{360})} \cdot \sin(2 \pi \cdot \frac{50}{360})
であるが、
 \sin(2 \pi \cdot \frac{50}{360})= \cos(2 \pi \cdot \frac{40}{360})
なのだから、これは 0 になる。 (証明終)

言い訳、その他

  • これは、「x が 30 度」だと見当がついていた場合の証明方法なので、角度が分からない場合だと使えないですね。


*1:三角形 ABC を 1/BC 倍すれば良い。

*2:偏角が負になることに注意。

数学の記号

数日前、友人から、
「数学の記号の『 \vdash 』ってどういう意味?」
と質問されました。
数学基礎論の授業でこの記号を見た覚えがあったので、「\vdash 基礎論」で google 検索をかけて、この記号が「証明可能」という意味だということを思い出したんだけど *1、これは結構ハードルが高い作業だと思います。


実際、「 \vdash 」の意味をネットで検索するためには、

  • 数式は TeX で書かれるのが一般的で
  •  \vdash」は TeX では「\vdash」と書くことで表示される

ということが分かっていないと出来ないはずです。
(誰か、これ以外のアイデアでこの記号の意味を調べることが出来た人はいますか?)


ぼくは前々から「数学と Web とは相性が悪いのではないか」と思ってるんだけど、「数学の記号や数式を検索するのが難しい (もしくは手間がかかる)」のも、その原因の一つな気がします。

追記 (1/22)

isao さんのコメント (2009/01/21 20:58) を受けて、数学の直和 *2 記号「 \oplus」でも似たことが出来ないか、調べてみました。
 \oplus」は「\oplus」と書くので (「O」に 「+」なので、そのまんまですよね!)、「\oplus 数学」で検索をかけても…、適当なのが出てこない!!
保健所の地図記号 *3 とほとんど同じなので、「数学記号 "保健所みたい"」で検索したら、見事に 0 件でした。。残念。


同様に、テンソル積の記号「\otimes」を「数学記号 "警察署みたい"」で google 検索しても無理でした。笑

*1:「\vdash 数学」でも行けますね。

*2:参考:直和 - Wikipedia

*3:参考:日本の地図記号の一覧 - Wikipedia

間違いがありました…。

先日書いた kmeans() 関数についての記事ですが、
kmeans 関数の 'singleton' オプション - ひらのの日記
間違いがありました。。すいません。
何も考えずにこう書いちゃったのですが↓

kmeans() のアルゴリズムは、朱鷺の杜↓
k-means法 - 機械学習の「朱鷺の杜Wiki」
に書かれているのと多分同じなので、以後、このサイトの表記を使うことにします。

実は Matlab だと、朱鷺の杜でいうところの「3」のアルゴリズムが微妙に違ってて、距離に "重み付け" がなされています。

誤:
全てのデータ  \mathbf{x}\in X を,各クラスタのセントロイド \mathbf{x}_i との距離  ||\mathbf{x}-\mathbf{x}_i|| を最小にするクラスタ X_i へ割り当てる


正: *1
データ  \mathbf{x}\in X を,(クラスタ X_i に対し定義される) 以下の評価関数  \varphi を最小にするクラスタX_i へ割り当てる.
 \varphi(X_i) := \alpha ||\mathbf{x}-\mathbf{x}_i||^2 (  x_i X_i のセントロイドとする. )
ここで, \alpha は以下のように定義される.

  1.  x が所属するクラスターが決まっていない場合: \alpha = 1
  2.  x が所属するクラスターが決まっている場合:まず,  x が所属しているクラスターを X_j とし, X_j の持つデータ数を  m とする. このとき, 以下のように定める.
    1.  i = j のとき: \alpha = m/(m - 1). ただし,  m = 1 のときは,  \alpha = 1.
    2.  i \neq j のとき: \alpha = m/(m + 1).

ソースコードを見たので、多分間違ってないと思うんだけど、、間違ってたら訂正します。


ぼくのような素人から見ると、ややこしいことやってるな〜と思うけど、こうやって、評価関数が、

設定されることで、所属するクラスターの変更が頻繁に起きるようになる、というメリットがあるのではないか…と個人的には思っています。


この 2 つの方法で結果が違ってくる例はいくらでも作ることが出来て、1 次元の場合でも、次のような例があります。

1 次元ユークリッド空間上に 4 点  p_1 = -2,  p_2 = 0,  p_3 = 0,  p_4 = 1 を取り、 X = \{p_1, p_2, p_3, p_4\} とする。
 X を、2 つのクラスタ X_1,  X_2 に分ける。
ここで、 X_1 の種  x_1 = -1 X_2 の種  x_2 = 1 とする。


このとき、クラスタリングの実行結果は以下のように異なってくる。

  • 「朱鷺の杜」版: X_1 = \{p_1, p_2, p_3\},  X_2 = \{p_4\}
  • Matlab 版: X_1 = \{p_1\},  X_2 = \{p_2, p_3, p_4\}



…それにしても、、こういうアルゴリズムって統一できないものなのでしょうかねぇ?
(もしくは、「違うことやってる!」ということがちゃんと分かる名前にするとか。)
「Kmeans」といったときに、人によっててんでばらばらのアルゴリズムを指している、というのはちょっとヤバイんじゃないか…と思います。
何となく、「宇治茶」騒動を連想してしまいます。
参考:宇治茶 - Wikipedia

*1:ユークリッド距離の場合です。

kmeans 関数の 'singleton' オプション

kmeans について何にも知らないのに、今回は kmeans について書いてみることにします。


Matlab では、kmeans 法を行う関数 kmeans() が準備されていて、
参考:Statistics and Machine Learning Toolbox Documentation
いろいろオプションも充実していて便利なようです。
kmeans() のアルゴリズムは、朱鷺の杜↓
k-means法 - 機械学習の「朱鷺の杜Wiki」
に書かれているのと多分同じなので、以後、このサイトの表記を使うことにします。


kmeans アルゴリズムを行うときは、セントロイドの初期値 *1 を決めなくてはいけないのですが、Matlab では、これをランダムにも手動ででも決めることが出来ます。
ただし、クラスターの種を変に設定すると、空のクラスターが出てくる場合があります。

クラスターが空になる例

1 次元ユークリッド直線上でクラスタリングを行う。
 X=\{0,1,3\}
を 2 つのクラスターに分類する。
クラスターの種を
 c_1=-1, c_2=0
とする。
このとき、《kmeans アルゴリズムの 3》 を行うと、
 X 上のどの点  p においても、
 ||p, c_1||>||p, c_2||
なので、
 X_1空集合になってしまう。


このままだと、《kmeans アルゴリズムの 2》で 0 で割り算をすることになってしまうので、通常はエラーとなってしまうのですが、「何とかしてクラスター数を減らさずに計算を続行したい!」というときのためのオプションが、タイトルにある 'singleton' オプションです。

kmeans([データ], [クラスター数], 'emptyaction', 'singleton' [ ,その他オプション])

とするだけです。
こうすると、空のクラスターがある空でないクラスターと入れ替わって、kmeans アルゴリズムが続行するのですが、恥ずかしながら最初は内部でどういう操作が行われているか全然意味が分かりませんでした。
マニュアルには

Action to take if a cluster loses all its member observations.
'singleton': Create a new cluster consisting of the one point furthest from its centroid.

と書かれてたんだけど、まず、『its centroid』がどれを指してるのか意味不明。
だから、『a new cluster consisting of the one point』がどれか良く分かっていませんでした。


『its centroid』ですが、ぼくは最初はてっきり、「空のクラスターのセントロイド」のことかと思ってました。
下図
(ここで、 c_i は、 X_i の種を表しています。初めて《kmeans アルゴリズムの 3》 を行った直後の状況です。)
で言うならば、『its centroid』とは、スッカラカンクラスタ X_1 の種  c_1 のことで、『そこから最も遠い点』というのは  p_6 のことかと思っていました。

だから、この操作で、 X_1 p_6 X_3 から奪って、

  •  X_1=\{ p_6 \}
  •  X_2=\{ p_1, p_2, p_3 \}
  •  X_3=\{ p_4, p_5 \}

となると思ってたんだけど、どうやらこういう意味ではないみたいです。
実際、この状況で、 c_2 とか  c_3 とかを動かして計算してみると、結果が変わってくるので、「空のクラスターのセントロイドと、各点との距離」という問題ではない、ということが分かります。


それで、いろいろ観察してみたんですけど、多分こういうアルゴリズムを行っているんじゃないかと思います。

  1. 任意の  X の点  p_i は、必ずどこかのクラスターに入っている。
  2. 従って、 p_i の入っているクラスターのセントロイドと、 p_i との距離は、各  p_i に対して一意に定まる。これを  d(p_i) と書く。
  3.  d(p_i) が最大となる  p_i が、求めるべき『the one point』である。複数個ある場合は、そのなかで添字の最も小さいものを選ぶ。

だから、上の図で言うならば、『the one point』は  p_3 であり、この操作で、 X_1 p_3 X_2 から奪って、

  •  X_1=\{ p_3 \}
  •  X_2=\{ p_1, p_2 \}
  •  X_3=\{ p_4, p_5, p_6 \}

となるわけです。


誤訳したぼくが悪いといえば全くその通りなんだけど、このマニュアルではやっぱり分かりにくいと思う。。

*1:以後、これを「クラスターの種」と呼ぶことにします。