最小値を与えるxの値は書くべきか?

†お気持ち表明†したくなったので書く

事前背景

例えば以下のような問題があったとする

$x$が実数のときの$y=x^2-4x$の最小値を求めよ

という問題があったとする。

このとき、最小値は-4である。ところで、この問題は「xがいつのときに最小値を取るか」について問われていない。しかしそのような状況であったとしても。最小値を取るxの値を答案に書くべきか?という問題である。

結論から言うと、もし書くとした場合、「なんのために書くべきか」をきちんと理解するべきである。

 

問題を抽象化すると以下の通りになる。

 

$A$を集合として、$f:A→\mathbb{R}$とする。このとき、$x\in A$における$f(x)$の最小値を求めよ

 

で、これを解くことを考える。そしてその答えが$m $であったとする。

 

答案というのは証明をするためのものなので、$f(x)$の最小値が$m $であることを示すためには以下の2つを示す必要がある。

  • 全ての$x\in A$について$f(x)\geq m $であること
  • ある$x_0\in A$が存在して$f(x_0)=m $であること

で、主題である「最小値を与える$x$の値は書くべきか?」についてだ・

 

「最小値を与える$x$の値を書く」というのは「ある$x_0\in A$が存在して$f(x_0)=m $であること」を示すための手段のうちの一つにすぎないのである。

 

冒頭の例に戻る

$x\in\mathbb{R}$における$f(x)=x^2-4x$の最小値を求める問題は、以下の2ステップで示すことになる。

1. $x^2-4x=(x-2)^2-4\geq -4$である。つまり任意の$x\in\mathbb{R}$について$f(x)\geq -4$である

2. $f(2)=-4$である。つまり、ある実数$x$が存在して$f(x)=-4$である。

 

2. についてだが、命題「ある実数$x$が存在して$f(x)=-4$」というのを示すために、そのような条件を満たす$x$を具体的に構成することで示したのである。「ある実数$x$が存在して$f(x)=-4$」という命題が出てきたとき、実際に条件を満たすものをドン!と出してくれば、それは紛れもなく真であることがわかる。

 

何度も言うが、「最小値を与える$x$の値を書く」というのは、「ある$x_0\in A$が存在して$f(x_0)=m $であること」を示すための手段のうちの一つである。

 

つまり、最小値を与える$x$の値を書かなかったとしても、実際に最小値を与えるような値が存在すること示すことができればそれでOKなのである。

 

そんな例あるのかよ?というわけで例題

$x\in\mathbb{R}$における、$f(x)=(e^x-3x)^2$の最小値を求めよ。

以下、最小値が0であることを示す。

1.については、$f(x)$は2乗の形をしているので、明らかに$f(x)\geq 0$が全ての実数$x$で成立することから成り立つ。

2.については、$f(x)=0$となるような実数$x$の存在をいえばOK

ここで、$g(x)=e^x-3x$とする。$f(x)=\{g(x)\}^2$が$x$について恒等的に成り立つことに注意

ここで、$g(x)$は任意の実数$x$で連続であり、$g(0)=1\gt 0,g(1)=e-3\lt 0$である。よって、中間値の定理より$0\lt c\lt 1$を満たして$g(c)=0$となるような実数$c$が存在する。このような$c$は実数であるため定義域に含まれており、さらに$f(c)=0$となる。つまり最小値を満たすような実数が存在することが示せたので、実際に0が最小値であることが示された。(証明終)

 

 

という感じで、$f(c)=0$の解が定義域上に存在することさえ示せてしまえば、実際にそのような値がどんな形をしているかは全く気にする必要はないのである。

 

高校数学においては、中間値の定理は証明されずfactとして与えられており、高校数学においては実際に解がどのように構成されるかということは全く扱わないのである。

このような状況では、最小値を与える$x$の値を書くことはなくなる。

つまり、このような状況においては「最小値を与える$x$の値を書く」以外の方法で「ある$x_0\in A$が存在して$f(x_0)=m $であること」を示せたので、これで問題ないということになる。

 

最後になるが、「最小値を与える$x$の値を書くべき」というのはあまり本質的ではなく、根底には「実際に最小値を達成するような$x$の値は存在するかどうかを明らかにする」というモチベが存在している。最小値を与える$x$の値を書くというのは、あくまでそのモチベを達成するための手段の一つというわけである。

斜回転体はパップスギュルダンで計算できる

斜回転体って難しいですよね。自分も大学入試の滑り止めを受けるときに出て返り討ちにされたので苦手意識があります。

斜回転体といえば傘型分割ですが、これはパップスギュルダンでもいけるのでは?と思ったので解いてみることにしました。

 

例題

問題は受験の月からパクってきました。

examist.jp

$y=x^2$と$y=x$で囲まれた部分の面積を$y=x$の周りに一回転してできる立体の体積は?

 

解答

$D=\{(x,y)\in\mathbb{R}^2;x^2\leq y\leq x,0\leq x\leq 1\}$

とする。$D$を$y=x$を軸に回転させたときの体積を求めればよい。

 

① $D$の面積を求める

$$\iint_{D} dxdy=\int_{0}^{1}\int_{x^2}^{x}dydx=\int_{0}^{1}x-x^2dx=\dfrac{1}{2}-\dfrac{1}{3}=\dfrac{1}{6}$$

 

② $D$の重心を求める

$$\iint_{D}x dxdy=\int_{0}^{1}\int_{x^2}^{x}xdydx=\int_{0}^{1}x^2-x^3dx=\dfrac{1}{3}-\dfrac{1}{4}=\dfrac{1}{12}$$

$$\iint_{D}y dxdy=\int_{0}^{1}\int_{x^2}^{x}ydydx=\int_{0}^{1}\left[\dfrac{1}{2}y^2\right]_{y=x^2}^{y=x}dx$$

$$=\dfrac{1}{2}\int_{0}^{1}x^2-x^4dx=\dfrac{1}{2}\left(\dfrac{1}{3}-\dfrac{1}{5}\right)=\dfrac{1}{15}$$

 

$D$の重心は

$$\left(\dfrac{\iint_D xdxdy}{\iint_D dxdy},\dfrac{\iint_D ydxdy}{\iint_D dxdy}\right)$$

なので、

$$\left(\dfrac{1/12}{1/6},\dfrac{1/15}{1/6}\right)$$

$$=\left(\dfrac{1}{2},\dfrac{2}{5}\right)$$

一応検算しておくと、$(1/2,2/5)$は,$1/4=1/2^2\leq 2/5\leq 1/2$なので$D$内部にあるので問題なさそうである。(一般的に重心が領域内部にあるという保証はないが、今回の例題の場合$D$は凸*1な領域なので、重心は内部にあると言えるはず。)

 

③重心と直線の距離

点と直線の距離の公式を使えばOK。点$(1/2,2/5)$と直線$x-y=0$との距離を求める。

$$\dfrac{|1\cdot 1/2 - 1\cdot 2/5+0|}{\sqrt{1^2+(-1)^2}}=\dfrac{\sqrt{2}}{20}$$

 

④体積

パップスギュルダンの定理は回転体の体積を求めるときに

(Dの面積)×2π×(Dの重心と直線の距離)

で求められるので、

答えは$$\dfrac{1}{6}\cdot 2\pi\cdot \dfrac{\sqrt{2}}{20}=\dfrac{\sqrt{2}\pi}{60}$$

 

感想

パップスギュルダンといえばドーナツの体積を計算して楽しんでるイメージしか無かったですが、こういう斜回転体にこそ真価を発揮するのではと思いました。

*1:凸の定義:任意のp_1,p_2 in D,t in [0,1]に対して、p_1 t+p_2 (1-t)がDの元となる

数独自動ソルバのしくみ

はじめに

数独というパズルがあります。まあこんな感じのものです。

 

数独の例(画像はwikipediaより) https://ja.wikipedia.org/wiki/%E3%83%95%E3%82%A1%E3%82%A4%E3%83%AB:Sudoku-by-L2G-20050714.svg


これをコンピュータに自動で解かせるにはどうすればいいでしょうか?

 

答えは†線形計画法†です。

アルゴリズム

$729=9^3$変数の線形計画法をします。具体的には

$$0\leq x_{i,j,k}\leq 1,(1\leq i,j,k\leq 9$$

です。これらの変数は、$(i,j)$マスが$k$なら1になって、そうでないならば0とします。

 

そして、1つのマスに1つの数字しかないので$\sum_{k} x_{i,j,k}\leq 1$です。すべてのi,jについて同じようなことをするので、条件式が$9^2=81$個あります

同じ行に同じ数字は2個以上ないので、$\sum_{j} x_{i,j,k}\leq 1$です。同上で条件式81個追加です。

同じ列に同じ数字は2個以上ないので、$\sum_{i} x_{i,j,k}\leq 1$です。同上で条件式81個追加です。

また、3x3のマスに同じ数字が2個以上ないので、

すべての$a,b\in\{1,4,7\},1\leq k\leq 9$について、$\sum_{i=0}^{2}\sum_{j=0}^{2}x_{i+a,j+b,k}\leq 9$

です。条件式81個追加

 

以上、条件式はすべて合わせると$81\times 4=324$個あります。

この条件のもとで、$\sum_{i,j,k}x_{i,j,k}$を最大化せよ。

という線形計画問題を解けばOKです。

 

追記

変数の範囲は、(i,j)マスがkであることが確定するならば$x_{i,j,k}=1$であり、特にそのような条件がないならば$0\leq x_{i,j,k}\leq 1$と設定します

 

pythonには線形計画問題を解いてくれるライブラリがあるので、それに丸投げします。

入力のフォーマットとしては、空欄の場所には代わりに0を入れてスペース区切りにすればOKです。

gist.github.com

 

inputの例

5 3 0 0 7 0 0 0 0
6 0 0 1 9 5 0 0 0
0 9 8 0 0 0 0 6 0
8 0 0 0 6 0 0 0 3
4 0 0 8 0 3 0 0 1
7 0 0 0 2 0 0 0 6
0 6 0 0 0 0 2 8 0
0 0 0 4 1 9 0 0 5
0 0 0 0 8 0 0 7 9

 

冒頭で紹介した数独の問題を投げると一瞬で答えを返してくれます。

 

プログラムの3行目をK=5にすれば25x25の数独にも対応してくれます。

ネットで適当に拾った問題を試しに入力したら7秒程度で答えが出力されました。

注意点

解が複数あるときにバグりそうです

余談

プログラムを大規模に書き換える必要がありますが、他の大抵の変わり種のナンバープレース線形計画法で殴れると思います。ただ立式が大変そうです

例えば、「$(i_1,j_1)$マスより$(i_2,j_2)$マスのほうが大きい」という条件を付け足したいなら

$$x_{i_1,j_1,k_1}+x_{i_2,j_2,k_2}\leq 0$$

を、すべての$(k_1\lt k_2)$に対して追加すればOKです。$N\times N$のナンプレならば条件を$\dfrac{N(N-1)}{2}$個追加することになります。

他にも、

「$(i_1,j_1)$マスの数字と$(i_2,j_2)$マスの数字は連続してる」

や、

「$(i_1,j_1)$マスの数字と$(i_2,j_2)$マスの数字は等しい」

や、

「$(i_1,j_1)$マスの数字と$(i_2,j_2)$マスの数字の和はxだ」

も同様にいくつかの条件を追加することで達成できると思います。

このプログラムにそこまでの価値があるのかはわかりませんが、勝手に使ったり改変したりして、どうぞ

懸賞問題でズルするのに使えそうですが、知りません。

最後に

この記事は「毎月1記事更新しよう」という目的をギリギリ達成させるために月末にやっつけで書かれています。

この記事ではscipyの線形計画法ライブラリを使っていますが、pulpを使っている先行研究を見つけたのでそれを見たほうがいいです。なんならpulpを使ったほうがコードを直感的に書けるので良さそうです。

地球上の2点における距離と角度の計算方法

概要

とあるプログラムを書くときに必要になったので。

設定

地球の半径は$R$

地球上の2点$P_i(i=1,2)$は緯度が$\theta_i$で経度が$\varphi_i$である。

このとき、

$$-\pi/2\lt \theta_i\lt \pi/2$$

$$-\pi\lt \varphi_i\lt \pi$$

とする。

原点$O$を地球の中心

z軸を地球の自転軸

つまり、$(0,0,R)$を北極点、$(0,0,-R)$を南極点として、

本初子午線と赤道の交点が$(R,0,0)$

赤道上・東経90°の地点が$(0,R,0)$

となるような座標系を考えてみる。

すると、

$$\overrightarrow{OP_i}=(x_i,y_i,z_i)=(R\cos{\theta_i}\cos{\varphi_i},R\cos{\theta_i}\sin{\varphi_i},R\sin{\theta_i})$$

距離

距離が近いなら、単にベクトルの差のノルムを求めたらそんなに誤差は出ないけど、2点の距離が離れるほど誤差がでかくなるので、きちんと実装しましょう。

地球のような球面だと、最短距離は大円を求めればOK

つまり$O,P_1,P_2$を通る平面で地球を切断したときの扇形$OP_1P_2$の孤の長さを求めればよい。

それは$\angle P_1OP_2$を求めればよいが、それは内積を計算すればいいことがわかる。

つまり、

$$R\cdot \cos^{-1}{\left(\dfrac{\overrightarrow{OP_1}\cdot \overrightarrow{OP_2}}{R^2}\right)}$$

が答え

 

この式は、P_1,P_2がどこにあっても正しく動くはず

角度

P_1から見たときのP_2の角度を求める。

「P_2から見たときのP_1の角度」ではないことに注意。後述

 

点$P_1=(x_1,y_1,z_1),x_1^2+y_1^2+z_1^2=R^2$に観測者がいるとき、

真北方向のベクトルは

$$\overrightarrow{N_1}=\left(-\dfrac{x_1z_1}{\sqrt{x_1^2+y_1^2}},-\dfrac{y_1z_1}{\sqrt{x_1^2+y_1^2}},\sqrt{x_1^2+y_1^2}\right)$$

真東方向のベクトルは

$$\overrightarrow{E_1}=\left(\dfrac{y_1R}{\sqrt{x_1^2+y_1^2}},-\dfrac{x_1R}{\sqrt{x_1^2+y_1^2}},0\right)$$

 

$\overrightarrow{E_1}$ベクトルと$\overrightarrow{N_1}$ベクトルは共にノルムが$R$であることに注意する。

(注:$P_1$が北極点または南極点である時、$\sqrt{x_1^2+y_1^2}=0$であるため分母に0が入って正常に計算できない。北極点から見たときはあらゆる方角が南、南極点から見たときはあらゆる方角が北であることに注意せよ)

このとき

$$Y=\overrightarrow{N_1}\cdot \overrightarrow{P_1P_2}$$

$$X=\overrightarrow{E_1}\cdot \overrightarrow{P_1P_2}$$

とすると、$(0,0)$から$(X,Y)$への方角と等しくなる

よって、ATAN2関数に$(X,Y)$を代入すれば答えが出る。

 

ここで、$P_2$がちょうど$P_1$の対蹠点にある場合については、$\overrightarrow{P_1P_2}$が$\overrightarrow{N_1},\overrightarrow{E_1}$と直交するので、$(X,Y)=(0,0)$となる。この調子でプログラムを組むと、ATAN2関数の仕様に沿った答えが出力される。ところで、この場合にはどの方向に進んでも$P_2$に到達できるので、当然ATAN2(0,0)の方向に進んでも目的地に到達できるわけである。よって、答えとしては間違ってはいないことになる。

 

結局、このプログラムを組んだ場合、$P_1$が北極点or南極点でない限り正しい答えが出る。

 

注意

「地点P_1から見た地点P_2の方角」

と、

「地点P_2から見た地点P_1の方角」

は必ずしも真逆になるとは限らない。

例えば東京から真東に進むと、途中で赤道を通過して南半球のチリに到達するのだが、途中の真東へ進む経路と赤道上の交点において、その地点から真西に進むと普通に赤道を一周してしまう。この場合、逆方向に進んだら東京には戻れない。



これ急に見て三角関数やんとはならないやろ

きっかけ

いつものようにボケーっとしながらツイッター見てたら面白そうなものを見つけた。

 

これ急に見て三角関数やんとはならないやろ

【画像】

 

画像の内容は、とある参考書を写したものである。

「$x^2+y^2=4$のときの$2x^2+3xy-y^2$の最大値・最小値を求めよ」という問題に対して、$x=2\cos{\theta},y=2\sin{\theta}$という変数変換をして解いているというような解法が書かれいた。

 

逆張りオタクをこじらせてしまったので、三角関数を使わない解法を考えてみることにした。ついでに一般化もしてみる。

考察

$n\in\mathbb{N}$として、$x=(x_1,\ldots,x_n)^\top \in\mathbb{R}^n$と書く。

とりあえず、$\parallel x\parallel =1$のときの$x^\top Ax$の最大値と最小値を求めよという問題を解くことにする。

ここで、$\parallel x\parallel:=\sqrt{x_1^2+\cdots+x_n^2}=\sqrt{x^\top x}$である。

また、$A$は$n$次元正方行列だが、

$$x^\top A x=\sum_{i,j}a_{i,j}x_i x_j$$

となっているが、

$$a_{i,j},a_{j,i}\to \dfrac{a_{i,j}+a_{j,i}}{2},\dfrac{a_{i,j}+a_{j,i}}{2}$$

という変換をしても式としては同じであるため、$A$は対称行列であるとしてもよい。

対称行列だと何が良いかというと、直交対角化ができる。

つまり、$A=P\Lambda P^{-1}$であり、

このとき$P$は直交行列、$\Lambda$は対角行列である。

直交行列なので、逆行列と転置が等しくなる。つまり、$P^\top=P^{-1}$となる。

また、対角行列$\Lambda$を成分表示すると以下のようになる。

$$\Lambda =\begin{pmatrix}\lambda_1&\cdots &0\\\vdots&\ddots&\vdots\\0&\cdots&\lambda_n\\\end{pmatrix}$$

ここで、$\lambda_1\geq \cdots\geq \lambda_n$としてもよい。(対称行列なので、固有値は全て実数であることが保証される。)

このとき、$y=P^{-1}x$と変数変換すると、

$1=x^{\top}x=x^{\top}PP^{\top}x=y^{\top}y$

である。よってこの問題は

$\parallel y\parallel =1$のときの$y^{\top}\Lambda y=\sum_{i=1}^{n}\lambda_i y_i^2$の最大値・最小値を求めれば良いことになる。

これは、

$$\lambda_n=\sum_{i=1}^{n}\lambda_n y_i^2\leq \sum_{i=1}^{n}\lambda_i y_i^2\leq \sum_{i=1}^{n}\lambda_1 y_i^2=\lambda_1$$

である。ここで、一個目の不等式は$y=(0,\ldots,0,1)^\top$のとき、二個目の不等式は$y=(1,0,\ldots,0)^\top$のときに等号が成立するので、実際に$\lambda_1,\lambda_n$がそれぞれ最大・最小となる。

ここで、$e_i$を$i$番目の基本ベクトルとする。$y=e_1$のときに最大となって$y=e_n$のときに最小となることに注意。このとき、$\Lambda e_i=\lambda_i e_i$なので、$\lambda_i, e_i$は$\Lambda $の固有値固有ベクトルの組となる。

ここで、$x^\top Ax$が最大・最小となるような$x$はそれぞれ$Pe_1,Pe_n$である。

ここで、

$$APe_i=PP^{-1}APe_i=P\Lambda e_i=\lambda_i Pe_i$$

なので、$\lambda_i,Pe_i$は$A$の固有値固有ベクトルの組となる。

つまり、$x^\top Ax$が実際に最大・最小となるような$x$は最大値・最小値(これらは$A$の固有値となる)に対応する固有ベクトルとなる。

 

これで終わりかと思いきや、件のツイートをよく見てみるとなんと、$x\geq 0,y\geq 0$という条件が課されているのである。

ここら辺の条件をなんとかするとしたら三角関数が一番手っ取り早いのかもしれないと思ったが、なんとかして逆張りを完遂したいので頑張る。

$n-1$次元球面上のある領域$D$について$\overline{D}:=D\cup \partial D$とする。

$x\in\overline{D}$における$x^\top Ax$の最大値・最小値を求めることにする。

ここで、$\overline{D}$はコンパクトなので、最大値・最小値を達成するような点は必ず存在する。最大値・最小値は無いというようなことは考えなくてもOK。

 

ここで、$\overline{D}$内に、$Pe_1$というベクトルが含まれているならば最大値は$\lambda_1$としても問題ないはずである。最小値も同じ

 

一方、含まれてない場合はどうすればいいかというと、これは見た感じ領域の境界点に最大最小を達成する点がありそうなので、これを示す。

ここで、今までは$x^\top Ax$の最大最小を考えていたが、以降では

$$f(x):=\dfrac{x^\top Ax}{x^\top x}$$

について考える。このとき、斉次化されているので$x$を定数倍しても値は変わらなくなっている。

($f(rx)=f(x),r\in(0,\infty) $)

よって、この場合、$D'=\{x\in\mathbb{R}^n\setminus \{0\}: \dfrac{x}{\parallel x\parallel}\in D\}$という領域内で考えてもいいことになっている。

ここで、

$\overline{D'}=\{x\in\mathbb{R}^n\setminus \{0\}: \dfrac{x}{\parallel x\parallel}\in \overline{D}\}$である。

ここで、この式を$x_i$番目の変数で偏微分すると

$$\dfrac{\partial}{\partial x_i}\dfrac{x^\top Ax}{x^\top x}=\dfrac{1}{(x^\top x)^2}\left( (x^\top x)\cdot \dfrac{\partial }{\partial x_i}(x^\top Ax)-(x^\top A x)\cdot \dfrac{\partial }{\partial x_i}(x^\top x)\right)$$

ここで、

$$\dfrac{\partial}{\partial x_i}(x^\top x)=2x_i$$

$$\dfrac{\partial}{\partial x_i}(x^\top A x)=\dfrac{\partial }{\partial x_i}(-a_ix_i^2+2\sum_{j=1}^{n}a_{i,j}x_ix_j)=2\sum_{j=1}^{n}a_{i,j}x_j=2(Ax)_i$$

となっているので、

$$\dfrac{\partial}{\partial x_i}f(x)=\dfrac{2}{(x^\top x)^2}\left( (x^\top x)\cdot (Ax)_i -(x^\top A x)\cdot x_i\right)$$

つまり、

$$\nabla f(x)=\dfrac{2(x^\top x)Ax-2(x^\top Ax)x}{(x^\top x)^2}$$

である。ここで、$y$が微小なベクトルであると仮定したとき、

$$f(x+y)\simeq f(x)+y^\top \nabla f(x)$$

である。ここで、$\nabla f(x)$が0ベクトルでないならば、ある微小なベクトル$y$を取れば$y^\top \nabla f(x)\gt 0$であり、$(-y)^\top \nabla f(x)\lt 0$となる。

このような状況において、$f(x-y)\lt f(x)\lt f(x+y)$となる。

もし$x\in D'$ならば、ある$\varepsilon\gt 0$が存在して$B(x,\varepsilon)\subset D'$なので、$x$は$f$を最適化しなくなる。つまり、ある微小なベクトル$y$が存在して、$f(x-y)\lt f(x)\lt f(x+y)$かつ、$x+y\in D',x-y\in D'$となるため、少なくとも$x$よりも大きく/小さくなる点が存在する。

 

では、$\nabla f(x)=0$となるのはどのような場合か?

それは、

$$(x^\top x)Ax=(x^\top Ax)x$$

であるようなときである。

これをよく見てみると、$Ax$と$x$が平行なベクトルであることが要求される。

つまり、$x$は$A$の固有ベクトルでなくてはいけない。

逆に$x$が固有ベクトルなら、固有値を$\lambda $とすると計算すれば$\nabla f(x)=0$が求まる。

結局、内部の点で最大値・最小値となるならば、それは固有ベクトルである。

 

$x$が固有ベクトルの場合

最大値・最小値を達成するような点が固有ベクトルであることがわかったが、それが最大固有値・最小固有値に対応する固有ベクトルかどうかはまだ言い切れない。

ここで、$x$は$A$の固有ベクトルとなるが、$x$に対応する固有値を$\lambda$とおく。ここで、$\lambda\neq \lambda_1,\lambda \neq \lambda_n$とする

このとき、

$$f(x+y)\simeq f(x)+y^\top H_f(x)y$$

となる。ここで、$H_f(x)$はヘッセ行列であり、$(i,j)$成分は$\dfrac{\partial}{\partial x_i}\dfrac{\partial}{\partial x_j}f(x)$である。

$$H_f(x)_{i,j}=\dfrac{\partial}{\partial x_j}\dfrac{\partial}{\partial x_i}\dfrac{x^\top Ax}{x^\top x}=\dfrac{\partial}{\partial x_j}\dfrac{1}{(x^\top x)^2}\left( 2(x^\top x)(Ax)_i-(x^\top A x)2x_i\right)$$

これを頑張って計算すると

$$=\dfrac{2}{(x^\top x)^3}\left(-2(x^\top x)((Ax)_i x_j+(Ax)_j x_i)+(x^\top x)^2 a_{i,j}+4(x^\top Ax)x_ix_j-(x^\top x)(x^\top Ax)\delta_{i,j}\right)$$

ここで、$\delta_{i,j}$はクロネッカーのデルタである($i=j$ならば$\delta_{i,j}=1$,$i\neq j$ならば$\delta_{i,j}=0$)

よって、ヘッセ行列は(注:$I_n$は単位行列

$$H_f(x)=\dfrac{2}{(x^\top x)^3}\left(-2(x^\top x)((Ax)x^\top +x(Ax)^\top)+(x^\top x)^2 A+4(x^\top Ax)xx^\top -(x^\top x)(x^\top Ax) I_n\right)$$

ここで、$Ax=\lambda x$であることから、

$$H_f(x)=\dfrac{2}{(x^\top x)}\left( A -\lambda I_n\right)$$

この行列の固有値は$\dfrac{2}{x^\top x}(\lambda_1-\lambda),\ldots,\dfrac{2}{x^\top x}(\lambda_n-\lambda)$である。

仮定より$\lambda_1-\lambda \gt 0,\lambda_n-\lambda\lt 0$

なので、$Pe_1$の方向に進めば、$f$の値はより大きくなり、$Pe_n$の方向に進めば、$f$の値はより小さくなる。よって、この場合は$x$において$f$は最大でも最小でもない。

よって、最大でない固有値に対応する固有ベクトルが最大値を達成することはない。最小についても同様。

 

結局、領域内部において最大値/最小値を達成するような点は最大固有値/最小固有値に対応する固有ベクトルであることがわかった。

で、もしそのような固有ベクトルが領域にない場合、その他の内部の点が最大値/最小値を達成することはない。よって、消去法により境界において最大・最小を達成することになる。

 

以上を踏まえると、最大/最小固有値に対応する固有ベクトルが領域内に存在しないならば、領域内における最大値/最小値については、境界の部分だけ探索すれば良いことになる。

結論

長々と書いたが、以下の結論が従う。

$x\in\mathbb{R}^n$として、$A$を$n$次対称行列として、$D\subset S^{n-1}$を領域とする。ここで、$S^{n-1}=\{x\in\mathbb{R}^n:\parallel x\parallel=1\}$である。

位相は$S^{n-1}$における部分位相を考えている。この位相において$D$は開集合であると仮定している。

このとき、$x\in \overline{D}$という範囲内での$x^\top Ax$の最大値・最小値について、以下が従う。($\overline{D}$はコンパクトなので最大値・最小値が必ず存在する。)

 

$A$の固有値は全て実数であることが保証されるが、その中で最大なものを$\lambda_1$とする。

もし、$\lambda_1$に対応する固有ベクトル$v_1$であって、$x\in\overline{D}$であるようなものが存在するならば、$v_1$において$x^\top Ax$は最大となる。

もしそのような固有ベクトルが存在しないならば、最大を達成するような$x$は領域の境界上のどれかである。

 

最小値についてもほぼ同じである。

 

$A$の固有値は全て実数であることが保証されるが、その中で最小なものを$\lambda_n$とする。

もし、$\lambda_n$に対応する固有ベクトル$v_n$であって、$x\in\overline{D}$であるようなものが存在するならば、$v_n$において$x^\top Ax$は最小となる。

もしそのような固有ベクトルが存在しないならば、最小を達成するような$x$は領域の境界上のどれかである。

 

 

解答

以上を踏まえて問題を解く

$$A=\begin{pmatrix}2&3/2\\3/2&-1\end{pmatrix}$$

として、

まずは$x^2+y^2=1$における最大・最小を求める。

 

これの固有方程式は$(2-t)(-1-t)-(3/2)^2=\dfrac{1}{4}(4t^2-4t-17)$なので、

これを解くことで$t=\dfrac{1\pm 3\sqrt{2}}{2}$が固有値であると分かる。

で、最大の固有値は$\dfrac{1+3\sqrt{2}}{2}$であり、最小の固有値は$\dfrac{1-3\sqrt{2}}{2}$である。

よって、$x\geq 0,y\geq 0$という条件を無視するならば、

$(x,y)=\pm \dfrac{1}{\sqrt{4+2\sqrt{2}}}(1+\sqrt{2},1)$のときに最大値$\dfrac{1+3\sqrt{2}}{2}$

$(x,y)=\pm \dfrac{1}{\sqrt{4+2\sqrt{2}}}(-1,1+\sqrt{2})$のときに最小値$\dfrac{1-3\sqrt{2}}{2}$

を達成するはずである。

 

ここで、$x\geq 0,y\geq 0$の条件を付け加えるとどうなるだろうか。

$(x,y)=\dfrac{1}{\sqrt{4+2\sqrt{2}}}(1+\sqrt{2},1)$とすると$x\geq 0,y\geq 0$を満たすので、このときに実際に最大値を達成する。よって最大値は$\dfrac{1+3\sqrt{2}}{2}$で問題ない。

しかし、最小値を達成するベクトルは$x\geq 0,y\geq 0$の範囲内に無い。

つまり、$x\geq 0,y\geq 0$という条件を追加した場合の最小値は領域の境界にあることになる。

ここで、境界は$(x,y)=(1,0),(0,1)$からなる2点集合となるので、それぞれについて計算すると、

$(x,y)=(1,0)$のとき$2$

$(x,y)=(0,1)$のとき$-1$

なので、最小値は$-1$である。($(x,y)=(0,1)$のときに最小値を達成する。)

 

以上が$x^2+y^2=1$のときの答えである。

ツイートにおける問題文は$x^2+y^2=4$となっている。これについては$(x,y)\to(2x,2y)$と変換すると$x^2+y^2=1$の場合と$x^2+y^2=4$の場合が一対一対応する。

よって$x^2+y^2=4$の場合を求めるには、$x^2+y^2=1$としたときの答えを単純に4倍すればよい。

よって、最大値は$2+6\sqrt{2}$であり、最小値は$-4$である。

 

おまけ(ラグランジュの未定乗数法)

$$f(t,x)=x^\top Ax-t x^\top x$$

というものを考える

$A$の$(i,j)$成分を$a_{i,j}$として、$A^\top=A$を仮定

ここで、

$$\dfrac{\partial}{\partial x_i}f(t,x)=\dfrac{\partial}{\partial x_i}\left(-a_{i,i}x_i^2+2\sum_{j=1}^{n}a_{i,j}x_ix_j -2tx_i^2\right)+(x_iに依存しない部分)$$

なので、

$$\dfrac{\partial }{\partial x_i}f(t,x)=2\sum_{j=1}^{n}a_{i,j}x_j-2tx_i$$

なので、

$$\nabla f(t,x)=2Ax-2tx=2(A-t I_n)x$$

となる。$\parallel x\parallel =1$で$\nabla f(t,x)=0$となるためには、

少なくとも$A-tI_n$が非正則行列でならなくてはいけない。このとき$t$は$A$の固有値となり、極値となるようなベクトル$x$は$A$の固有値$t$に対する固有ベクトルとなる。つまり、最大・最小を達成するようなベクトルも固有ベクトルとなる。あとは同じ

 

あとがき

ツイートの問題を解くだけなら、極座標変換したほうが一番手っ取り早い気がした。

でも、多次元だと極座標変換したら面倒くさくなるし対称性も崩れるのでこの記事でやったように線形代数で殴ったほうが良いのかもしれない。

 

余談

行列のノルムの定義(のうちの1つ)は以下の通りである。

$$\parallel A\parallel :=\sup_{\parallel x\parallel=1}\parallel Ax\parallel$$

ここで、

$$\parallel Ax\parallel=\sqrt{x^\top A^\top Ax}$$

このとき、$A^\top A$は対称行列であることに注意

つまり、$\parallel x\parallel =1$の場合の$x^\top A^\top Ax$の最大値を求めればよいことになる。

よって、行列のノルムは$A^\top A$の最大の固有値のルートを取ったものとなる。

ここで、任意の$x\in\mathbb{R}^n$について$x^\top A^\top Ax=\parallel Ax\parallel ^2\geq 0$であるため、$A^\top A$は半正定値行列となる。よって全ての固有値が非負であることが保証されるので安心してルートを取ることができる。

F_q係数多項式をq乗すると…?

はじめに

任意の$\mathbb{F}_2$係数多項式$f$において

$$f(x^2)=\{f(x)\}^2$$

が成り立つ。これは比較的簡単に示せる。

$$f(x)=a_0+\cdots+a_{n-1}x^{n-1}$$

なら、

$$\{f(x)\}^2=\sum_{i=0}^{n-1}a_i^2x^{2i}+2\sum_{i\lt j} a_ia_jx^{i+j}$$

となるが、$\mathbb{F}_2$の標数は2なので、第二項は消えるのである。また、$a_i\in\mathbb{F}_2$ならば$a_i^2=a_i$となるため、恒等式が成立する。

 

これを見るとある疑問が浮かぶ。一般化できないだろうか?と。

つまり、

$\mathbb{F}_q$係数多項式$f$は常に$f(x^q)=\{f(x)\}^q$となるだろうか?

結論から言うとこれは正しいのだが、上記の2乗の場合の方法を一般化するとなると、ゴリゴリ計算することになって面倒くさいので、もっといい方法でやってみることにする。

事前知識

有限体について

$\mathbb{F}_q$を位数$q$の有限体とする。しかし、有限体の位数は必ず素数の整数乗となるため、$q=p^m $とする。ここで、$p$は素数であって、$m $は正の整数である。

このとき、$\mathbb{F}_q$は標数$p$、つまり$p$倍すると消えることに注意である。

また、任意の$x\in\mathbb{F}_q$に対して、$x^q=x$となる。これは有名事実である。

ちなみに、$\mathbb{F}_q$は$\mathbb{F}_p$の有限次拡大で拡大次数は$m $である。

二項係数modについて

$q=p^m $,$p$は素数,$m $は正の整数とする。

このとき、$1\leq k\leq q-1$に対して、

$${}_qC_{k}\equiv 0\pmod{p}$$

となる。

【証明】

証明の方法は色々ある。一番まともな方法は、ルジャンドルの定理を使う方法である。

ルジャンドルの定理については下記リンク参照

manabitimes.jp

${}_qC_{k}=\dfrac{q!}{k!(q-k)!}$なので、${}_qC_{k}$が$p$で割り切れる回数は、

$$\sum_{n=1}^{\infty}\lfloor\dfrac{p^m}{p^n}\rfloor-\lfloor\dfrac{k}{p^n}\rfloor-\lfloor\dfrac{p^m-k}{p^n}\rfloor$$

となる。

ここで、ガウス記号の性質として、

$$\lfloor x+y\rfloor \geq \lfloor x\rfloor +\lfloor y\rfloor$$

というものがある。三角不等式と不等号の向きが逆なことに注意

よって、この無限級数(実際は有限項)の各項はすべて非負となる。

ここで、$n=m $のとき、

$$\lfloor\dfrac{p^m}{p^m}\rfloor-\lfloor\dfrac{k}{p^m}\rfloor-\lfloor\dfrac{p^m-k}{p^m}\rfloor=1-0-0=1\gt 0$$

なので、この級数の値は正の整数となる。よって$p$で割り切れる回数が0でないため、$p$の倍数である。

 

【別解】

他にも、飛び道具を使う方法もある。Lucasの定理というものを思い出す。

Lucasの定理についてはリンク参照

manabitimes.jp

${}_qC_{k}$を$p$で割った余りを考えるとき、$q,k$をp進数で表記ものとする。

ここで、

$$q=a_0+a_1p+\cdots +a_mq^m,k=b_0+b_1p+\cdots+b_mq^m $$

である。ただし、$0\leq a_i,b_i\lt p$である。

このとき、

${}_qC_{k}\equiv \prod_{i=0}^{m}{}_{a_i}C_{b_i}\pmod{p}$

である。しかし、$q=p^m $なので、$a_0=a_1=\cdots=a_{m - 1}=0,a_m=1$である。

また、$k\lt p^m $より$b_m=0$である。

よって、${}_{a_m}C_{b_m}={}_{1}C_{0}=1$に注意すると

${}_qC_{k}\equiv \prod_{i=0}^{m - 1}{}_{0}C_{b_i}\pmod{p}$

となる。ここで、$0\lt k$より、ある$i$について$b_i\neq 0$となる。

このとき${}_{0}C_{b_i}\equiv 0\pmod{p}$となる。よって、$m $個の因子の中に0が必ずあるのでmod pで0になる。

証明

証明は帰納法でやる。以下の場合で証明できたら帰納法がきちんと走るはずである。

  1. $f(x)=a$(定数)の場合
  2. $f(x)=x$の場合
  3. $f(x),g(x)$で成り立つならば、$f(x)g(x)$でも成り立つ
  4. $f(x),g(x)$で成り立つならば、$f(x)+g(x)$でも成り立つ

2.と3.で$f(x)=g(x)=x$とすれば$x^2$でも成り立つ。

同様に、$x^{n-1}$で成り立つならば$f(x)=x^{n-1},g(x)=x$とすれば$x^n$でも成立する。よって、任意の自然数$n$について$f(x)=x^n$は条件を満たす。

$f(x)=x^n,g(x)=a$で3.を適用すれば$ax^n$でも成り立つので、すべての単項式で条件を満たす。

任意の$n-1$次以下の多項式で成り立つならば、

$f(x)=a_0+\cdots+a_{n-1}x^{n-1},g(x)=a_nx^n$はそれぞれ条件をみたすので、

$f(x)+g(x)=a_0+\cdots+a_{n-1}x^{n-1}+a_nx^n$でも条件をみたすので、任意の$n$次以下の多項式で成り立つ。

(また、1.より、任意の0次以下の多項式=定数で成り立つ)

よって、任意の多項式で条件を満たす。

 

それぞれを示す。

1.$f(x)=a$(定数多項式)のとき

$$f(x^q)=a,\{f(x)\}^q=a^q$$

となる。ここで、$a\in\mathbb{F}_q$なので、$a^q=a$である。よって

$$a=f(x^q)=\{f(x)\}^q=a^q$$

となる、

 

2. $f(x)=x$のとき

$$f(x^q)=x^q,\{f(x)\}^q=x^q$$

よって明らかに成り立つ。

 

3. $f(x),g(x)$で成立$\Rightarrow$ $f(x)g(x)$で成立

$$f(x^q)g(x^q)=\{f(x)\}^q\{g(x)\}^q=\{f(x)g(x)\}^q$$

つまり、

$$(fg)(x^q)=\{f(x)\}^q\{g(x)\}^q=\{(fg)(x)\}^q$$

より成立

 

4. $f(x),g(x)$で成立$\Rightarrow$ $f(x)+g(x)$で成立

二項定理より、

$$\{f(x)+g(x)\}^q=\sum_{i=0}^{q}{}_{q}C_{i} \{f(x)\}^i \{g(x)\}^{q-i}$$

となる。しかし、さきほど示した定理より、${}_{q}C_{i}$は$0\lt i\lt q$ならば$p$の倍数である。$\mathbb{F}_q$の標数は$p$であるため$p$の倍数なら0とみなしても良い。

よって、

$$\{f(x)+g(x)\}^q=\sum_{i=0}^{q}{}_{q}C_{i} \{f(x)\}^i \{g(x)\}^{q-i}\equiv \{f(x)\}^q+\{g(x)\}^q=f(x^q)+g(x^q)$$

よって、

$$\{(f+g)(x)\}^q=(f+g)(x^q)$$

となる。

 

1.2.3.4すべて示せたので、帰納法が無事走って、任意の多項式$f$で$f(x^q)=\{f(x)\}^q$となることが示された。

おまけ(多変数多項式

実は、任意の$\mathbb{F}_q$係数多変数多項式$f$で

$$f(x_1^q,\ldots,x_k^q)=\{f(x_1,\ldots,x_k)\}^q$$となる。

証明は1変数の場合とほぼ同じ。より詳しく説明すると

  1. $f(x_1,\ldots,x_k)=a$(定数)の場合
  2. $f(x_1,\ldots,x_k)=x_i$の場合$(i=1,\ldots,k)$
  3. $f(x_1,\ldots,x_k),g(x_1,\ldots,x_k)$で成り立つならば、$(fg)(x_1,\ldots,x_k)$でも成り立つ
  4. $f(x_1,\ldots,x_k),g(x_1,\ldots,x_k)$で成り立つならば、$(f+g)(x_1,\ldots,x_k)$でも成り立つ

を示せばOKである。

1.2.3より単項式で成り立つことがわかって、4.を加えることで任意の多項式でも成り立つことが分かる。

おまけ2(有理式の場合)

有理式の場合でも成り立つことを示す。

有理式は多項式の商体なので、

任意の多項式$f(x_1,\ldots,x_k),g(x_1,\ldots,x_k)$について、$(f/g)(x_1,\ldots,x_k):=f(x_1,\ldots,x_k)/g(x_1,\ldots,x_k)$で成り立つことが示せればOKである。

ここで、

$$(f/g)(x_1^q,\ldots,x_k^q)=\dfrac{f(x_1^q,\ldots,x_k^q)}{g(x_1^q,\ldots,x_k^q)}$$

$$=\dfrac{\{f(x_1,\ldots,x_k)\}^q}{\{g(x_1,\ldots,x_k)\}^q}$$

$$=\left(\dfrac{f(x_1,\ldots,x_k)}{g(x_1,\ldots,x_k)}\right)^q$$

$$=\{(f/g)(x_1,\ldots,x_k)\}^q$$

より成り立つ。

 

二項係数 mod 素数のあまり見かけない計算方法

atcoder.jp

 

これを頑張って計算すると以下の式(を998244353で割った余り)を出力すればいいことがわかります。

$$\sum_{M=0}^{N}\binom{N}{M}\cdot \binom{K+M+1}{N+1}$$

この式の導出方法はこの記事の本筋からそれるので省略します。解説参照!

さて、制約によると$K\leq 10^{18},N\leq 3\times 10^5$でmodの法は998244353(素数)みたいです。これを時間内に求めるにはどうすればいいでしょう???

 

うまくいかない方法1

二項係数のmod 素数を求める典型テクニックは、

$1!,2!,...,n!$と、$1/1!,...,1/n!$のテーブルを$O(n)$で作ってmod逆元を組み合わせる方法が一般的である。しかしこの方法の場合、$K+N+1$くらいの大きさのテーブルを作る必要があるが、そんなクソデカテーブルなど作れるわけがないのでこの方法ではうまくいきません。

 

うまくいかない方法2

$$\binom{K+M+2}{N+1}=\dfrac{K+M+2}{K+M+1-N}\binom{K+M+1}{N+1}$$

となるため、$\binom{K+1}{N+1}$を求めてから順々に$\binom{K+M+1}{N+1}$を計算するとうまくいきそうに見えるかもしれない。

しかし、この場合、$K$が998244353よりも大きくなる場合があるので、$K+M+1-N$が998244353の倍数となる可能性がある。そうなるとゼロ除算になってしまうため、この方法だと正常に計算できない。よって、うまくいかない。

 

うまくいかない方法3

judge.yosupo.jp

$m $が大きすぎて間に合いません!!!

結局、Lucasの定理を使おうにも、998244353個分の階乗modを記録する必要があり、それは無理である。

うまくいかない方法4

要は$\binom{K+M+1}{N+1}$を出力するクエリを高速化すればよい。

で、これはつまり

$$(K+M+1)\cdot (K+M)\cdot \cdots \cdot (K+M-N+1)$$

を求めれればよい。

つまり、

$$ a_i=(K - N + 1)\cdot (K - N + 1 + 1)\cdot \cdots \cdot (K - N + i) $$

という順列を前計算で取得して、

$$\dfrac{a_{N+M+1}}{a_{M}}$$

を計算すればよいのでは?

と、思いたくなるが、この場合も、$K-N+1,\ldots,K+N+1$の中に998244353の倍数が混じっていた場合正しく計算できなくなる。よって、うまくいかない。

 

うまくいく方法の例

実は、うまくいかない方法4は実はいい線をいっている。

区間積の方法を累積積で計算したのがまずかった。以下の方法だとうまくいく。

長さ$2N+1$のセグメントツリーを定義する。ここで、

・初期値は$K-N+1,\ldots,K+N+1$という長さ$2N+1$の列

演算子は$x\otimes y=(x*y) \mod 998244353$

単位元は(もちろん)1

とする。

このセグメントツリーに対して、$M $ごとに$a_M\otimes \cdots \otimes a_{N+M}$を計算して出力する。

 

するとクエリ一回あたりの計算量は$O(\log{N})$となり、最終的に全体の計算量は$O(N\log{N})$となるので通る。

また、$K-N+1,\ldots,K+N+1$の間に998244353の場合があっても正しく動く。

セグメントツリーは任意のモノイドに対して動く。

モノイドとは、結合法則を満たして単位元が存在すればよかった。

つまり、逆元が存在しなくても問題ないので、区間内に998244353の倍数が存在しても正しく動くのである。