凸二次優(yōu)化問題

Theory. 設(shè)

是實對稱正定矩陣，，則求解凸二次優(yōu)化問題

等價于求解線性方程組

。

Proof. 二次型二階可導(dǎo)，極小值點處梯度為零，現(xiàn)對優(yōu)化的目標(biāo)函數(shù)求梯度。二次型本質(zhì)上具有：

計算梯度的分量表達(dá)式：

合在一起寫成矩陣形式：

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

最速梯度下降法

稱優(yōu)化目標(biāo)函數(shù)的梯度為殘量（Residue），即是當(dāng)前解對于線性方程組的不滿足量：

由于函數(shù)是凸函數(shù)，極值點一定存在。當(dāng)前解處函數(shù)的梯度值表示了函數(shù)值上升最快的方向（梯度方向上方向?qū)?shù)最大）。那么沿著相反的方向每迭代一步就會更加靠近最優(yōu)的極小值解。于是，我們定義迭代關(guān)系：

其中

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

在最速下降法中，將

作為自變量代入原函數(shù)，并看作的函數(shù)，對其進(jìn)行最小化：

同樣的，駐點處梯度為零。根據(jù)鏈?zhǔn)角髮?dǎo)法則：

最后一步由于

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

算例設(shè)計與實驗

參考南科大的數(shù)值分析作業(yè)題EH計算設(shè)計SUSTech的算例一（由于版權(quán)關(guān)系無傳送門）

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

：

其中

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

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

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

推導(dǎo)條件數(shù)計算公式 矩陣

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

那么其特征值介于

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

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

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

，x0表示解的初始向量。偽隨機數(shù)算法參數(shù)給定如下面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;

這里設(shè)定最優(yōu)解的分量在

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

這里設(shè)定算法的停機標(biāo)準(zhǔn)為當(dāng)前梯度變?yōu)榻咏诹愕臉O小量：

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

實驗結(jié)果

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

算法變體

這里清川對原始算法做了一些小改進(jìn)，在

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

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

這是因為這里選取的歐式范數(shù)（

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

總結(jié)

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

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。