采用梯度下降優化器(Gradient Descent optimizer)結合禁忌搜索(Tabu Search)求解矩陣的全部特征值和特征向量

[前言]

---

　　對於矩陣(Matrix)的特征值(Eigens)求解，采用數值分析(Number Analysis)的方法有一些，我熟知的是針對實對稱矩陣(Real Symmetric Matrix)的特征值和特征向量(Characteristic Vectors)求解算法——雅克比算法（Jacobi）。Jacobi算法的原理和實現可以參考[https://blog.csdn.net/zhouxuguang236/article/details/40212143]。通過Jacobi算法可以以任意精度近似求解任意實對稱矩陣的所有特征值和特征向量。Jacobi算法利用了實對稱矩陣的特性得到了迭代求解公式，而對於一般方陣，不能利用Jacobi求解特征值和特征向量。因此更一般的求解方法是基於矩陣的特征值和特征向量的定義計算得到的。對於尚不清楚或者已經遺忘矩陣的特征值和特征向量的定義的讀者，建議參考知乎[https://www.zhihu.com/question/21874816]，每一個回答都從某個角度講述了矩陣的特征根和特征向量的意義。需要注意的是，學習矩陣的性質和變換，不得不深入其內在的幾何意義。同時，矩陣和群論(Galois Group Theory)之間也可以建立起聯系：在求解矩陣特征值時，要求行列式的值為0。因為，所以，這樣問題就轉化為一個一元（該元為）高次多項式求零點的問題。而關於這個方程是否有根式解的判定，需要借助Galois群論知識，主要是通過置換轉化為可解群的過程。當然，如果該矩陣不是實對稱矩陣，那么該一元高次方程有可能在復數域存在解，意味着該矩陣可能存在復數的特征值和特征向量，其求解過程將會變得稍微復雜。

　　本文的研究問題是：如何采用近似最優化算法求解得到全局的所有可能解。

　　本文的實驗對象是：一般方陣在實數域內的特征值和特征根的求解算法。

　　本文的研究思路是：對梯度下降算法進行改進，通過引入禁忌表、局部最優解的感知和跳出機制實現持續解空間搜索，進而得到全部全局最優解。我們順帶研究了如何避免梯度下降過慢的解決方案。

[實驗]

---

 　　采用梯度下降方法一般只能求得一個近似最優解，這無法滿足我們尋求所有滿足的和的要求，因此，我們需要對梯度下降方法進行優化，使之能迭代求得所有滿足的解。另外，如果我們的方法幸運地成功了，那么意味着，我們不但能夠擴展到一般多解問題的求解上，同時還能對全局最優解的近似求解方法進行改進。這樣想的話，我們的工作將變得十分有趣！

　　那么接下來，我們采用Python語言在Tensorflow這個十分好用的機器學習框架下搭建我們的機器學習模型，並利用Tensorflow中已為我們集成的tf.train.GradientDescentOptimizer這個梯度下降優化器進行目標求解。代碼如下：

 

# python2.7+, tensorflow 1.5rc+
import numpy as np
import numpy.linalg as lin
import tensorflow as tf

''' Only REAL NUMBER is considered as a solution,
    solutions in complex field is nonsense.
    For a real symmetric matrix, all the eigen 
    values are real, and all the eigen vectors
    are real. Every pair of the eigen vector
    are orthogonal.
'''

''' Solving eigen values of a matrix using
    numpy linear algorithm tools.
'''
def NP_Eigen(mat):
    assert mat.shape[0]==mat.shape[1]
    return lin.eig(mat)

''' Normalizing vector with L2 metrics. 
'''
def TF_L2_Norm(m):
    t_scale = tf.sqrt(tf.reduce_sum(m*m))
    t_sign = m[0]/tf.abs(m[0])
    t_scale = 1.0/tf.maximum(t_scale, 1e-7)
    tf.stop_gradient(t_scale)
    tf.stop_gradient(t_sign)
    return t_sign*t_scale*m

''' Solving eigen values of a matrix using
    stochastic gradient descent optimizer.
    mat: the matrix to solve
    N: the maximum number of iteration
    M: the maximum number of jump
    R: the effective radius of tabu spots
'''
def GD_Eigen(mat, N, M, R):
    assert mat.shape[0]==mat.shape[1]
    eigens = np.zeros([mat.shape[0]])
    vectors = np.zeros(mat.shape)
    t_A = tf.constant(mat)
    t_lambda = tf.Variable(tf.random_uniform([1])-0.5)
    t_x = tf.Variable(tf.random_uniform([mat.shape[1],1])-0.5)
    t_x = TF_L2_Norm(t_x)
    t_err_pure = tf.matmul(t_A, t_x) - (t_lambda*t_x)
    t_err_pure = tf.reduce_mean(t_err_pure*t_err_pure)
    t_reg = tf.constant(0, tf.float32)
    opt = tf.train.GradientDescentOptimizer(learning_rate=0.1) 
    sess = tf.Session()
    sess.run(tf.global_variables_initializer())  
    for counter in xrange(mat.shape[0]):
        if counter>0:
            t_dist = tf.reduce_mean(tf.square(tf.reshape(t_x,[t_x.shape[0]*t_x.shape[1]])-vectors[counter-1]))
            t_cond = tf.cast(t_dist<R*R, tf.float32)
　　　　　　　t_reg_add = 1.0/tf.maximum(t_dist, 1e-7)
            t_reg_add = t_cond*t_reg_add
            t_reg = t_reg + t_reg_add
        t_err = t_err_pure + t_reg
        trainer = opt.minimize(t_err)
        err = 1e7
        err_last = 0
        i = 0
        j = 0
        while i<N and j<M:
            i += 1
            err_last = err
            _, err, reg = sess.run([trainer, t_err, t_reg])
            if err < 1e-12:
                print 'Root#', counter, ' found at iter#', i
                print 'Solution is: ', sess.run([t_lambda, t_x])
                break
            elif abs(err_last - err) < 1e-15:
                print 'Trapped! Reinitialize vars and search again!'
                sess.run(tf.global_variables_initializer())
                i = 0
                j += 1
            elif err/(err_last-err)>N-i:
                print 'Optimize too slow at speed: ', err_last-err
                sess.run(tf.global_variables_initializer())
                i = 0
                j += 1
            #print err, reg
        # saving this solution
        eigens[counter], x = sess.run([t_lambda, t_x])
        vectors[counter] = x.reshape([mat.shape[1]])
    return eigens, vectors.T

def main():
    tf.set_random_seed(128)
    A = np.array([[1,3,-4],[3,-2,5],[-7,6,0]], dtype=np.float32)
    A = A.dot(A.T)
    A = A/np.sqrt(np.sum(A*A))
    s, u = GD_Eigen(A, 10000, 100, 1e-1)
    print s
    print u
    s, u = NP_Eigen(A)
    print s
    print u

main()

　　注意到，在上述代碼中，我們采用了numpy.linalg(Numpy Linear Algorithm Library)這一個工具庫中的eigen函數對同一目標矩陣求解了特征值和特征向量，eigen函數采用的精確解的求解方法，因此其結果總是正確的，為此，我們將我們的代碼運行得到的結果和eigen函數的結果對比發現，我們得到的所有特征值和特征向量都是正確的，雖然他們都或多或少存在1e-5精度上的誤差，但這不影響我們的正確率。另外，即便你更改上述程序中的初始隨機種子(tf.set_random_seed(128))，依然能夠觀測到類似結果。這表明我們的算法是正確的。

　　那么，接下來我將慢慢解釋上述代碼的算法由來。

　　1. 梯度下降算法和一些近似最優求解算法一樣，無法一次性求得全部滿足最優化目標的解，而且，有時候得到的解並非全局最優，而僅僅是局部最優。其原因是梯度下降算法只是朝着當前最佳的誤差下降的方向搜索解空間。除了精確解算法，在非凸目標函數的最優化過程中，貪婪規則下的算法很容易就陷入到局部最優解，而且接下來會困死在那里，因為局部最優解所在的空間附近是平坦的(在這一點上對誤差的梯度的模趨近0)，而且稍微動一下就會出現上坡路(在該點的鄰域內任意點，其對誤差的梯度方向背離該點，這意味着誤差要想下降，必須沿着梯度反方向——也就是朝着該點的方向運動，這使得該點成為了一個局部穩定點，小擾動根本不會使之逃離該點)，這使得Agent就仿佛陷入到了一個陷阱里，按照貪婪規則再也爬不出來了。 為此，我們提出局部最優解的跳出機制，但在這之前，我們需要一個方法能監測Agent是否陷入到了局部最優解中。其次我們才能針對這種情形設計迫使Agent跳出的機制。

　　2. 為了搜索得到全部滿足優化目標的解，我們需要設計持續搜索的機制，並且為了避免陷入循環或浪費搜索時間，我們應該設計避免Agent回到已搜索過的最優解上。

　　3. 有時候，Agent的出發點不一樣會導致陷入不同的局部最優解，也可能會導致梯度下降的坡度太緩，優化速度過慢。因此，需要設計檢測下降速度的機制並針對該情形設計跳出機制。

　　4. 對於優化目標，我們需要設計正規化(Normalization)方法使得誤差的值域是合理的。而且這一點對於包含正則項的優化目標來說尤為重要，因為它決定正則系數的選取。

　　5. 由於持續搜索機制，我們需要考慮如何避免由之而來的新問題：新增部分對優化目標的最終解是否有影響？影響多大？能否完全避免？

 

　　幸運的是，上述的代碼中，這些問題我們已經全部解答了，而且在實驗中已經證實上述問題都得到了不錯的解決。

　　可以先看下運行日志：

 

Root# 0  found at iter# 1076
Solution is:  [array([0.05244394], dtype=float32), array([[0.6913288 ],
       [0.6959271 ],
       [0.19429319]], dtype=float32)]
Trapped! Reinitialize vars and search again!
Optimize too slow at speed:  6.117625e-08
Trapped! Reinitialize vars and search again!
Root# 1  found at iter# 403
Solution is:  [array([0.31607816], dtype=float32), array([[ 0.6778219 ],
       [-0.53152364],
       [-0.5079764 ]], dtype=float32)]
Trapped! Reinitialize vars and search again!
Optimize too slow at speed:  3.1590462e-06
Optimize too slow at speed:  6.239861e-08
Optimize too slow at speed:  6.187474e-08
Optimize too slow at speed:  3.874302e-06
Optimize too slow at speed:  6.367918e-08
Optimize too slow at speed:  6.7055225e-08
Optimize too slow at speed:  4.425645e-06
Optimize too slow at speed:  3.3862889e-06
Optimize too slow at speed:  3.0882657e-06
Optimize too slow at speed:  6.030314e-08
Optimize too slow at speed:  6.1816536e-08
Optimize too slow at speed:  1.356937e-06
Optimize too slow at speed:  6.088521e-08
Optimize too slow at speed:  6.251503e-08
Optimize too slow at speed:  9.071082e-07
Optimize too slow at speed:  6.6822395e-08
Optimize too slow at speed:  6.0070306e-08
Optimize too slow at speed:  4.5895576e-06
Optimize too slow at speed:  3.571622e-07
Optimize too slow at speed:  9.313226e-10
Optimize too slow at speed:  3.3304095e-06
Optimize too slow at speed:  6.589107e-08
Optimize too slow at speed:  6.0498714e-06
Optimize too slow at speed:  6.164191e-08
Optimize too slow at speed:  5.9138983e-08
Optimize too slow at speed:  6.0594175e-08
Optimize too slow at speed:  6.0827006e-08
Root# 2  found at iter# 382
Solution is:  [array([0.94728214], dtype=float32), array([[ 0.25024486],
       [-0.48287624],
       [ 0.839171  ]], dtype=float32)]
[0.05244394 0.31607816 0.94728214]
[[ 0.69132882  0.67782187  0.25024486]
 [ 0.69592708 -0.53152364 -0.48287624]
 [ 0.19429319 -0.50797641  0.83917099]]
[0.9472825  0.05244393 0.3160783 ]
[[-0.2502432   0.69133323 -0.6778176 ]
 [ 0.48287496  0.6959237   0.53152794]
 [-0.8391723   0.19428991  0.50797766]]

 

　　從上述運行日志中可以看到，Agent隨着持續搜索進程的進行，發現的解越多，以后搜索到新解的難度越大，這一方面是因為優化目標本身的解空間中：容易發現的解具有更強大的吸引子，它能把絕大多數隨機出發點上的Agent吸引到該解上。而不易被發現的解的吸引子的引力場較弱，大部分隨機出發點上的Agent難以找尋到通往該解的最快的路徑。

　　在我們設計的算法中，每找到一個解，我們都會將該解添加到禁忌表中，並且將該解作為禁區，設置一個反引力子迫使任何進入該禁區的Agent因為誤差梯度突然變得極大而被排擠出。該引力子的添加其實就是在優化目標上添加一項新的正則項，該正則項在有限范圍內（R的設置，在本實驗中取0.1）給Agent一個很大的懲罰。這樣Agent將不會停留在該禁區中，從而實現了持續搜索以及跳出機制。