ざっくりとMajorization techniqueを解説(参考: Modern Multidimensional Scaling)

##1. Majorization(優関数法)とは
　Majorizationとは関数(損失関数など)の直接的な最適化ができない際に，以下の条件を満たす補助関数に置き換え，それを繰り返し最適化することで元の関数を最適化する手法です．多変量解析の手法では，多次元尺度構成法(Multidimensional Scaling, MDS)や多次元展開法(Multidimensional unfolding)などの目的関数を最小化する際に利用されています．(発想はEMアルゴリズムにちょっと似てる?)

補助関数(auxiliary function, あまり厳密ではないです)
元の関数を$f(x)$とし，この記事内では損失関数とする．補助関数の性質は以下のようなものである．
1.元の関数$f(x)$より最小化を行うの簡単な形である．
2. 補助関数$g(x,y)$は，元の関数$f(x)$と比較して常に大きいか，もしくは等しい．つまり$f(x) \leq g(x,y)$を常に満たす
3. supporting pointと呼ばれる点をyとすると，補助関数はsupporting pointで元の関数$f(x)$に触れる．つまり，supporting point yで，$f(y) = g(y,y)$である

　性質2から，補助関数は_Majorization function_とも呼ばれます．

##2. どうして_Majorization function_を最適化するの?
　なぜ，_Majorization function_を最適化すると元の関数が最適化されるかという疑問は，性質2,3から導かれる以下の不等式を確認することで解決します．

$$f(x^* ) \le g(x^*, y) \le g(y, y) = f(y)$$

ここで，$x^* $は点yにおける補助関数$g(x,y)$をxに関して最小化する点を表しています．この不等式(sandwich inequalityと呼ばれます)から，両サイドの元の関数$f(x)$のみに注目すると，$f(x^* ) \le f(y)$が残り，補助関数を最適化することで元の関数の値が減少していくことがわかります．

##3. 具体的な例
　例えばMDSの損失関数(3行目は2行目の各項を置き換えただけです)

\begin{eqnarray}
stress(X) &=& \sum_{i<j}w_{ij}(\delta_{ij}-d_{ij}(X))^2\\
&=&\sum_{i<j}w_{ij}\delta_{ij}^2 + \sum_{i<j}w_{ij}d_{ij}^2(X)-2\sum_{i<j}w_{ij}\delta_{ij}d_{ij}(X)\\
&=&\eta_{\delta}^2 + \eta^2(X)-2\rho(X)
\end{eqnarray}

を最適化する問題を考えます．_stress_を行列形式で書くと，

\begin{eqnarray}
stress &=& \eta_{\delta}^2 + trX^TVX - 2tr\ X^TB(X)X
\end{eqnarray}

ここで$B(X)$の$(i,j)$成分は

\begin{eqnarray}
 b_{ij} &=& 
\left\{
\begin{array}{c}
&- \frac{w_{ij}\delta_{ij}}{d{ij}(X)}\quad &\mbox{for  i≠j and }d_{ij}(X) \neq0\\
 &0\quad  &\mbox{for i≠j and }d_{ij}(X)=0
\end{array}
\right.
\\b_{ii} &=& - \sum_{j=1,j\neq i}^n b_{ij}
\end{eqnarray}

としています．stress_において$\eta{\delta}^2$は定数，第2項は2次，第3項はちょっと複雑です．これをXに関して最適化する際に_Majorization function_として

\tau(X,Y) = \eta_{\delta}^2 + trX^TVX - 2tr\ X^TB(Y)Y

用います．この_Majorization function_は_stress_と比較して，第3項が線形となっており，全体としては2次の関数となっています．よって，この関数を最小化するX(反復計算におけるXの更新式)は

\nabla\tau(X,Y) = 2VX - 2B(Y)Y

を0とおくことで，

X_{update} = V^{+}B(Y)Y

と求めることができます．ここで，$V^+$は$V$のムーアペンローズ逆行列を表しています．$Y$に前のステップでの$X$の推定値を用いて_stress_が収束するまで繰り返し更新していくことでXの最適値を計算することができます．

　以下にMDSをスクラッチで実装したコードを載せます(蛇足)．

import numpy as np 
import matplotlib.pyplot as plt 

def mds(Dissimilarities,Weights,ndim = 2,limitIter=1000,nstart=15):
	'''
	Modern Multidimensional Scaling p.187-192

	Dissimilarities : symmetric, nonnegative and hollow(=all diagonal elements = 0)
					  I denote Dissimilarities as Del, abbreviation of Delta
	Weights : if dij, i-th and j-th elements of Dissimilarities, is missing, coordinate wij = 0.
			  otherwise w = 1
	ndim : the dimension of configuration
	'''
	Del = Dissimilarities.astype("float64"); W = Weights.astype("float64"); n = Del.shape[0]
	Xu = np.zeros((n,ndim)); ones = np.ones(n); J = np.identity(n) - (1/n)*np.outer(ones,ones.T)
	k = 0; 

	#calculate eta_delta V
	eta_del = 0; V = np.zeros((n,n))
	for i in range(n):
		for j in range(n):
			if i >= j:
				continue
			ei = np.zeros(n); ei[i] = 1
			ej = np.zeros(n); ej[j] = 1
			Aij = np.outer((ei-ej),(ei-ej).T)
			V += W[i,j] * Aij
			eta_del += W[i,j] * (Del[i,j]**2)
	print(eta_del); #print("V = ",V) #for check
	Vminv = np.linalg.inv(V + np.outer(ones,ones.T)) - (1/(n**2))*np.outer(ones,ones.T)

	'''
	We MUST repeat the implemantations starting from below, 
	because we use random start configuration. 
	'''

	#initial configration
	X_initial = np.random.randn(n*ndim).reshape(n,ndim);X = np.copy(X_initial);
	Z = np.copy(X)
	print("initial configuration\n",X_initial)
	Bx = np.zeros((n,n))
	for i in range(n):
		for j in range(n):
			dx = np.sqrt(np.sum((X[i,]-X[j,])**2))
			if dx == 0:
				Bx[i,j] = 0
			else:
				Bx[i,j] = -(W[i,j]*Del[i,j])/dx
	for i in range(n):
		Bx[i,i] = -np.sum(np.hstack([Bx[i,:i],Bx[i,(i+1):]]))
	#print("B(X)",Bx) #for check
	stress = eta_del + np.trace(np.dot(X.T,np.dot(V,X))) - 2*np.trace(np.dot(X.T,np.dot(Bx,X)))
	stress_old = stress; diff = stress_old - stress
	stress_list = []; stress_list.append(stress)

	while k == 0 or (diff > 10**(-6) and k <= limitIter):
		k += 1
		print(k," : ",stress)
		

		#Z = np.copy(X); 
		Bz = np.zeros((n,n))
		for i in range(n):
			for j in range(n):
				dz = np.sqrt(np.sum((Z[i,]-Z[j,])**2))
				if dz == 0:
					Bz[i,j] = 0
				else:
					Bz[i,j] = -(W[i,j]*Del[i,j])/dz
		for i in range(n):
			Bz[i,i] = -np.sum(np.hstack([Bz[i,:i],Bz[i,(i+1):]]))
		#print("B(Z)",Bz)
		Xu = np.dot(Vminv,np.dot(Bz,Z))

		Bxu = np.zeros((n,n));
		for i in range(n):
			for j in range(n):
				dxu = np.sqrt(sum((Xu[i,]-Xu[j,])**2))
				if dxu == 0:
					Bxu[i,j] = 0
				else:
					Bxu[i,j] = -(W[i,j]*Del[i,j])/dxu
		for i in range(n):
			Bxu[i,i] = -np.sum(np.hstack([Bxu[i,:i],Bxu[i,(i+1):]]))

		stress = eta_del + np.trace(np.dot(Xu.T,np.dot(V,Xu))) - 2*np.trace(np.dot(Xu.T,np.dot(Bxu,Xu)))
		stress_list.append(stress)
		diff = stress_old - stress

		if diff < 0:
			print("!!diff is negative!!")
			break
		Z = Xu
		stress_old = stress

	#print(Xu) #for check
	return Xu, stress_list, X_initial


if __name__ == "__main__":
	Disimi = np.array([[0,11,2,3,4,6],[11,0,3,5,9,7],[2,3,0,8,2.5,1],
		[3,5,8,0,1.7,5],[4,9,2.5,1.7,0,12],[6,7,1,5,12,0]])

	W = np.ones((6,6))

	Xhat, list_s, X_initial =  mds(Disimi,W,ndim=2)
	print(" Given dissimilarities\n",Disimi)

	# Visualization 
	names = ["a","b","c","d","e"]
	# values of loss function
	plt.plot(np.arange(1,(len(list_s)+1)),list_s,linewidth=4,
		color="green",linestyle="dashdot",marker="D")
	plt.title("Decreesing loss values");
	plt.xlabel("number of iteration");plt.ylabel("loss value")
	plt.show()

	# initial configration
	for (i,j,name) in zip(X_initial[:,0],X_initial[:,1],names):
		plt.plot(i,j,'*')
		plt.annotate(name,xy=(i,j))
	plt.title("initial configration")
	plt.show()

	# estimated configration
	for (i,j,name) in zip(Xhat[:,0],Xhat[:,1],names):
		plt.plot(i,j,'D')
		plt.annotate(name,xy=(i,j))
	plt.title("estimated configration")
	plt.show()

(アルゴリズム通り実装しており，重複する処理等を削ったりしていないので無駄が多いコードですが，問題なく動きます．)
##4. 参考書籍
　Borg and Groenen(2005). Modern Multidimensional Scaling. Springer.を参考にしました．
　MDSは気が向いたら，より丁寧な記事をアップロードします．