凸二次優化問題

Theory. 設

是實對稱正定矩陣，，則求解凸二次優化問題

等價于求解線性方程組

。

Proof. 二次型二階可導，極小值點處梯度為零，現對優化的目標函數求梯度。二次型本質上具有：

計算梯度的分量表達式：

合在一起寫成矩陣形式：

顯然，凸二次優化問題的極值條件等價于該線性方程組。凸二次優化問題在建模中十分常見，這說明討論線性方程組的求解方法具有普遍的實用價值。然而對于規模較大的問題，使用線性代數中的克萊姆法則暴力展開將導致時間開銷巨大，而高斯消元法算法流程又較為復雜。本文將介紹一種常見的數值分析方法，求得線性方程組的數值近似解。

最速梯度下降法

稱優化目標函數的梯度為殘量（Residue），即是當前解對于線性方程組的不滿足量：

由于函數是凸函數，極值點一定存在。當前解處函數的梯度值表示了函數值上升最快的方向（梯度方向上方向導數最大）。那么沿著相反的方向每迭代一步就會更加靠近最優的極小值解。于是，我們定義迭代關系：

其中

表示迭代第輪的解，表示每輪迭代的步長，即每一步下降多少的權重。之所以需要設計一個與相關的步長，是因為隨著迭代的進行，梯度是變化的。越靠近駐點的梯度將會越小，直到接近于0時收斂。同時還要考慮在接近于駐點時應該放緩步伐，否則會出現梯度在正負之間頻繁震蕩，解在最優解左右搖擺的情況。

在最速下降法中，將

作為自變量代入原函數，并看作的函數，對其進行最小化：

同樣的，駐點處梯度為零。根據鏈式求導法則：

最后一步由于

是正定矩陣，所以分母也是一個正定二次型，值不為零，可以直接除過來。由此，我們最終得到了解的迭代關系：

算例設計與實驗

參考南科大的數值分析作業題EH計算設計SUSTech的算例一（由于版權關系無傳送門）

生成實對稱正定矩陣 采用文獻[1]中的類似方法生成矩陣

：

其中

和分別為Householder矩陣和對角矩陣：

Householder矩陣是對稱正交矩陣，這時，對角矩陣的特征值就是

的特征值，參考矩陣的正交分解。取：

推導條件數計算公式 矩陣

的最大特征值和最小特征值分別為：

那么其特征值介于

之間，對稱正定矩陣的條件數公式為

于是可以反解處條件數的計算公式：

為了能夠重復試驗數據，使用線性同余法計算偽隨機數進行向量的初始化（也可以使用相同的隨機種子，調用編程語言內部的隨機值庫，這里尊重南科大的原題要求）。令xopt表示最優解

，x0表示解的初始向量。偽隨機數算法參數給定如下面MATLAB程序所示：v = zeros(n, 1); xopt = zeros(n, 1); x0 = zeros(n, 1); t = 0; for j = 1: n t = mod(t * 31416 + 13846, 46261); v(j) = t * (2 / 46261) - 1; end; t = 0; for j = 1: n t = mod(t * 42108 + 13846, 46273); xopt(j) = t * (5 / 46273) + 5; end;for j = 1: n t = mod(t * 42108 + 13846, 46273); x0(j) = t * (5 / 46273) - 10; end;

這里設定最優解的分量在

之間，不為零的初始向量分量在之間。根據以上公式可以根據向量得到矩陣的值，其后使用計算得到向量的值。這時將和作為問題的給定參數，x0作為初始向量，使用梯度下降的方式計算xk，并計算它與最優解的真值之間的相對誤差：

這里設定算法的停機標準為當前梯度變為接近于零的極小量：

Python程序代碼

import numpy as npdef initialize(cond, numb):gamma = (np.cos(np.pi / (numb + 1)) + cond * np.cos(numb * np.pi / (numb + 1)) - (cond - 1)) / (cond - 1)diags = np.array(list(map(lambda i: np.cos(i * np.pi / (numb + 1)) + 1 + gamma, np.arange(1, numb + 1))))SIGMA = np.diag(diags)v = np.zeros((numb, 1))xopt = np.zeros((numb, 1))x0 = np.zeros((numb, 1))t = 0for j in range(numb):t = (t * 31416 + 13846) % 46261v[j] = t * (2 / 46261) - 1t = 0for j in range(numb):t = (t * 42108 + 13846) % 46273xopt[j] = t * (5 / 46273) + 5for j in range(numb):t = (t * 42108 + 13846) % 46273x0[j] = t * (5 / 46273) - 10V = np.identity(numb) - 2 * v @ v.T / np.linalg.norm(v) ** 2A = V @ SIGMA @ V.Tb = A @ xoptreturn A, b, xopt, x0def gradientDescent(A, b, x0, epsilon):x_now = x0; epoch = 0gnorm = ginit = np.linalg.norm(A @ x0 - b)while gnorm / ginit > epsilon:g_now = A @ x_now - bgnorm = np.linalg.norm(g_now)xnext = x_now - (gnorm ** 2) / (g_now.T @ A @ g_now) * g_now # * LR_linearDecay(epoch)x_now = xnextepoch += 1return epoch, x_nowdef LR_linearDecay(epoch):return - epoch / 5e4 + 0.9def LR_expertDecay(epoch):return np.exp(- np.log(0.9) * (epoch / 4e3 - 1))if __name__ == '__main__':print("{:>8}".format("EPOCH/ERR"), end=" ")for n in range(1, 6):print("n={:>9}".format(100*n), end=" ")print()for c in range(3, 7):print(f"COND=1e+{c}", end=" ")for n in range(1, 6):A, b, xopt, x0 = initialize(cond=10**c, numb=100*n)epoch, xk = gradientDescent(A, b, x0, epsilon=1e-6)error = np.linalg.norm(xk - xopt) / np.linalg.norm(x0 - xopt)print("{:>5}/{:.4f}".format(epoch, error), end=" ")print()

MATLAB程序代碼

Linear_init.m

function [A, b, xopt, x0] = linear_init(cond, numb)gamma = (cos(pi / (numb + 1)) ...+ cond * cos(numb * pi / (numb + 1)) - (cond - 1)) / (cond - 1);diags = zeros(numb, 1);for i = 1: numb;diags(i) = cos(i * pi / (numb + 1)) + 1 + gamma;end;SIGMA = diag(diags);v = zeros(numb, 1); xopt = zeros(numb, 1); x0 = zeros(numb, 1);t1 = 0; t2 = 0;for j = 1: numbt1 = mod(t1 * 31416 + 13846, 46261); v(j) = t1 * (2 / 46261) - 1; end; for j = 1: numbt2 = mod(t2 * 42108 + 13846, 46273); xopt(j) = t2 * (5 / 46273) + 5; end;for j = 1: numbt2 = mod(t2 * 42108 + 13846, 46273); x0(j) = t2 * (5 / 46273) - 10; end;V = eye(numb) - 2 * (v * v') / norm(v, 2)^2;A = V * SIGMA * V';b = A * xopt; end

Linear_solu.m

function [epoch, x_now] = linear_solu(A, b, x0, epsilon)x_now = x0; epoch = 0;gnorm = norm(A * x0 - b, 2);ginit = gnorm;while gnorm / ginit > epsilong_now = A * x_now - b;gnorm = norm(g_now, 2);xnext = x_now - (gnorm^2) / (g_now' * A * g_now) * g_now;x_now = xnext;epoch = epoch + 1;end; end

Linear_impr.m

function [epoch, x_now] = linear_impr(A, b, x0, epsilon)lr = @(epoch) - epoch / 5e4 + 0.9;x_now = x0; epoch = 0;gnorm = norm(A * x0 - b, 2);ginit = gnorm;while gnorm / ginit > epsilong_now = A * x_now - b;gnorm = norm(g_now, 2);xnext = x_now - lr(epoch) * (gnorm^2) / (g_now' * A * g_now) * g_now;x_now = xnext;epoch = epoch + 1;end; end

Linear_eval.m

function linear_eval(type) % linear_eval('normal') for normal gradient descent algorithm % linear_eval('improv') for improved radient descent algorithmepsilon = 1e-6;fprintf('%8s ', 'EPOCH/ERR');for n = 1: 5fprintf(' n=%3d ', 100 * n);end;fprintf('n');for c = 3: 6fprintf('COND=1e+%d ', c);for n = 1: 5[A, b, xopt, x0] = linear_init(10^c, 100 * n);if strcmp(type, 'improv')[epoch, xk] = linear_impr(A, b, x0, epsilon);else[epoch, xk] = linear_solu(A, b, x0, epsilon);end;error = norm(xk - xopt, 2) / norm(x0 - xopt, 2);fprintf('%5d/%.4f ' , epoch, error);end;fprintf('n');end; end

實驗結果

1 EPOCH/ERR n= 100 n= 200 n= 300 n= 400 n= 500 COND=1e+3 3382/0.0786 8910/0.0682 14400/0.0511 15750/0.0492 15696/0.0525 COND=1e+4 3382/0.0786 8910/0.0682 14400/0.0511 15750/0.0492 15696/0.0525 COND=1e+5 3382/0.0786 8910/0.0682 14400/0.0511 15750/0.0492 15696/0.0525 COND=1e+6 3382/0.0786 8910/0.0682 14400/0.0511 15750/0.0492 15696/0.0525 -epoch / 5e4 + 0.9 EPOCH/ERR n= 100 n= 200 n= 300 n= 400 n= 500 COND=1e+3 203/0.0786 459/0.0682 419/0.0510 588/0.0494 544/0.0536 COND=1e+4 221/0.0786 534/0.0682 855/0.0513 538/0.0494 761/0.0535 COND=1e+5 203/0.0786 333/0.0682 617/0.0509 538/0.0494 761/0.0535 COND=1e+6 203/0.0786 534/0.0682 617/0.0509 588/0.0494 544/0.0536np.exp(- np.log(0.9) * (epoch / 4e3 - 1)) EPOCH/ERR n= 100 n= 200 n= 300 n= 400 n= 500 COND=1e+3 287/0.0786 462/0.0682 513/0.0512 579/0.0489 545/0.0535 COND=1e+4 313/0.0786 531/0.0682 705/0.0512 617/0.0496 710/0.0528 COND=1e+5 287/0.0786 412/0.0682 619/0.0510 617/0.0496 710/0.0528 COND=1e+6 287/0.0786 531/0.0682 619/0.0510 579/0.0489 545/0.0535

算法變體

這里清川對原始算法做了一些小改進，在

前面加上了一個人為設定的二次權重值。實驗結果中每部分的第一行代表了的取值方式，表中斜線前面代表停機時的總迭代次數，后面代表相對誤差。可以看出，在保持解的精度不變的前提下，加上人工修正的權重使得算法快了30倍。一開始清川想改進成隨機梯度下降，但顯然他理解錯了，SGD算法的隨機指的是神經網絡訓練時每輪迭代選取隨機的樣本，和這里沒有什么關系。但是清川從SGD中采納了指數衰減和線性衰減的方法，并加以調參就得到了改進。

這個改進并不能讓人喜悅，顯然它是在當前數據（矩陣A與b的值）下過擬合的。對于新的數據，該算法未必能夠加速。但是這說明了另一個問題：最速下降法并不是最速的。

這是因為這里選取的歐式范數（

分子上的二范數）并不一定是最合適的衡量標準。每輪迭代總是假設梯度在當前迭代時是近似不變的，以初始位置的梯度模擬整個步長上的梯度，這樣并不準確。而清川加了人工調參的二次權重后更好地模擬了梯度的變化，所以更快。如果想進一步實驗這一點，可以將直接用替換，去調參，預計仍然能得到一個很好的收斂速度。做實驗前最好對進行歸一化，否則調參過程中可能出現數值溢出。關于最速下降法不總是最速的這個結論，也可以參考梯度下降法和最速下降法區別。注意本文并未對這兩個概念作區分。

總結

以上是生活随笔為你收集整理的qr分解求线性方程组_梯度下降求解线性方程组算例设计的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。