優化三對角線係數矩陣的A * x = B解

Question

我有一個A*x = B形式的方程組，其中[A]是一個三對角線係數矩陣。使用Numpy求解器numpy.linalg.solve我可以求解x的方程組。

請參閱下面的示例以瞭解如何開發三角形[A] martix。該{B}載體，解決了x：

# Solve system of equations with a tridiagonal coefficient matrix 
# uses numpy.linalg.solve 

# use Python 3 print function 
from __future__ import print_function 
from __future__ import division 

# modules 
import numpy as np 
import time 

ti = time.clock() 

#---- Build [A] array and {B} column vector 

m = 1000 # size of array, make this 8000 to see time benefits 

A = np.zeros((m, m))  # pre-allocate [A] array 
B = np.zeros((m, 1))  # pre-allocate {B} column vector 

A[0, 0] = 1 
A[0, 1] = 2 
B[0, 0] = 1 

for i in range(1, m-1): 
    A[i, i-1] = 7 # node-1 
    A[i, i] = 8  # node 
    A[i, i+1] = 9 # node+1 
    B[i, 0] = 2 

A[m-1, m-2] = 3 
A[m-1, m-1] = 4 
B[m-1, 0] = 3 

print('A \n', A) 
print('B \n', B) 

#---- Solve using numpy.linalg.solve 

x = np.linalg.solve(A, B)  # solve A*x = B for x 

print('x \n', x) 

#---- Elapsed time for each approach 

print('NUMPY time', time.clock()-ti, 'seconds') 

所以我的問題涉及到上面的例子中的兩個部分：

1. 因爲我處理了[A]一個角矩陣，也稱爲帶狀矩陣，有沒有更有效的方法來解決方程組而不是使用numpy.linalg.solve？
2. 此外，有沒有更好的方法來創建三角對角矩陣，而不是使用for-loop？

上述示例根據time.clock()函數在Linux上運行約0.08 seconds。

的numpy.linalg.solve功能工作正常，但我試圖找到一種方法，在進一步加快解決方案的希望採取的[A]的三對角形式的優勢，然後將該方法更加複雜的例子。

Answer 1

有兩種立即性能改進（1）不使用循環，（2）使用scipy.linalg.solve_banded()。

我會寫代碼的東西更像

import scipy.linalg as la 

# Create arrays and set values 
ab = np.zeros((3,m)) 
b = 2*ones(m) 
ab[0] = 9 
ab[1] = 8 
ab[2] = 7 

# Fix end points 
ab[0,1] = 2 
ab[1,0] = 1 
ab[1,-1] = 4 
ab[2,-2] = 3 
b[0] = 1 
b[-1] = 3 

return la.solve_banded ((1,1),ab,b) 

可能有更優雅的方式來構建矩陣，但這個工程。

在ipython中使用%timeit原始碼花了112毫秒爲m = 1000。這個代碼在m = 10,000時需要2.94 ms，這是一個數量級更大的問題，但仍然快了兩個數量級！我沒有耐心等待m = 10,000的原始代碼。原來的大部分時間可能在構建陣列，我沒有測試這個。無論如何，對於大型數組來說，只存儲矩陣的非零值會更有效率。

Answer 2

有一個scipy.sparse矩陣類型，稱爲scipy.sparse.dia_matrix，它很好地捕獲了矩陣的結構（它將存儲3個數組，在「位置」0（對角線），1（上方）和-1（下方））。使用這種類型的矩陣，你可以嘗試scipy.sparse.linalg.lsqr來解決。如果你的問題有一個確切的解決方案，它會被發現，否則它會找到最小二乘意義的解決方案。

from scipy import sparse 
A_sparse = sparse.dia_matrix(A) 
ret_values = sparse.linalg.lsqr(A_sparse, C) 
x = ret_values[0] 

然而，這可能不是在剝削triadiagonal結構方面完全優化，有可能使這一快的理論方法。這種轉換對你的作用是將矩陣乘法開銷減少到必要條件：僅使用3個頻段。這與迭代求解器lsqr的組合應該已經產生加速。

注意：我並不是主張scipy.sparse.linalg.spsolve，因爲你的矩陣轉換爲csr格式。但是，用spsolve代替lsqr值得一試，尤其是因爲spsolve可以綁定UMFPACK，請參閱相關doc on spsolve。另外，請看this stackoverflow question and answer relating to UMFPACK

Answer 3

您可以使用scipy.linalg.solveh_banded。

編輯：你爲你的矩陣不是對稱的，我認爲這是CAN NOT使用上述。然而，如上面的評論中提到的，Thomas算法是爲這個偉大的

a =  [7] * (m - 2) + [3] 
b = [1] + [8] * (m - 2) + [4] 
c = [2] + [9] * (m - 2) 
d = [1] + [2] * (m - 2) + [3] 

# This is taken directly from the Wikipedia page also cited above 
# this overwrites b and d 
def TDMASolve(a, b, c, d): 
    n = len(d) # n is the numbers of rows, a and c has length n-1 
    for i in xrange(n-1): 
     d[i+1] -= 1. * d[i] * a[i]/b[i] 
     b[i+1] -= 1. * c[i] * a[i]/b[i] 
    for i in reversed(xrange(n-1)): 
     d[i] -= d[i+1] * c[i]/b[i+1] 
    return [d[i]/b[i] for i in xrange(n)] 

此代碼不能優化也不會使用np，但如果我（或任何這裏的其他精細人的）有時間，我會編輯它，以便它做到這些。目前，m = 10000時約爲10毫秒。

Answer 4

這可能會幫助 有一個函數creates_tridiagonal將創建三對角矩陣。還有另一個函數可以按照SciPy solve_banded函數的要求將矩陣轉換爲對角有序形式。

import numpy as np  

def lu_decomp3(a): 
    """ 
    c,d,e = lu_decomp3(a). 
    LU decomposition of tridiagonal matrix a = [c\d\e]. On output 
    {c},{d} and {e} are the diagonals of the decomposed matrix a. 
    """ 
    n = np.diagonal(a).size 
    assert(np.all(a.shape ==(n,n))) # check if square matrix 

    d = np.copy(np.diagonal(a)) # without copy (assignment destination is read-only) error is raised 
    e = np.copy(np.diagonal(a, 1)) 
    c = np.copy(np.diagonal(a, -1)) 

    for k in range(1,n): 
     lam = c[k-1]/d[k-1] 
     d[k] = d[k] - lam*e[k-1] 
     c[k-1] = lam 
    return c,d,e 

def lu_solve3(c,d,e,b): 
    """ 
    x = lu_solve(c,d,e,b). 
    Solves [c\d\e]{x} = {b}, where {c}, {d} and {e} are the 
    vectors returned from lu_decomp3. 
    """ 
    n = len(d) 
    y = np.zeros_like(b) 

    y[0] = b[0] 
    for k in range(1,n): 
     y[k] = b[k] - c[k-1]*y[k-1] 

    x = np.zeros_like(b) 
    x[n-1] = y[n-1]/d[n-1] # there is no x[n] out of range 
    for k in range(n-2,-1,-1): 
     x[k] = (y[k] - e[k]*x[k+1])/d[k] 
    return x 

from scipy.sparse import diags 
def create_tridiagonal(size = 4): 
    diag = np.random.randn(size)*100 
    diag_pos1 = np.random.randn(size-1)*10 
    diag_neg1 = np.random.randn(size-1)*10 

    a = diags([diag_neg1, diag, diag_pos1], offsets=[-1, 0, 1],shape=(size,size)).todense() 
    return a 

a = create_tridiagonal(4) 
b = np.random.randn(4)*10 

print('matrix a is\n = {} \n\n and vector b is \n {}'.format(a, b)) 

c, d, e = lu_decomp3(a) 
x = lu_solve3(c, d, e, b) 

print("x from our function is {}".format(x)) 

print("check is answer correct ({})".format(np.allclose(np.dot(a, x), b))) 


## Test Scipy 
from scipy.linalg import solve_banded 

def diagonal_form(a, upper = 1, lower= 1): 
    """ 
    a is a numpy square matrix 
    this function converts a square matrix to diagonal ordered form 
    returned matrix in ab shape which can be used directly for scipy.linalg.solve_banded 
    """ 
    n = a.shape[1] 
    assert(np.all(a.shape ==(n,n))) 

    ab = np.zeros((2*n-1, n)) 

    for i in range(n): 
     ab[i,(n-1)-i:] = np.diagonal(a,(n-1)-i) 

    for i in range(n-1): 
     ab[(2*n-2)-i,:i+1] = np.diagonal(a,i-(n-1)) 


    mid_row_inx = int(ab.shape[0]/2) 
    upper_rows = [mid_row_inx - i for i in range(1, upper+1)] 
    upper_rows.reverse() 
    upper_rows.append(mid_row_inx) 
    lower_rows = [mid_row_inx + i for i in range(1, lower+1)] 
    keep_rows = upper_rows+lower_rows 
    ab = ab[keep_rows,:] 


    return ab 

ab = diagonal_form(a, upper=1, lower=1) # for tridiagonal matrix upper and lower = 1 

x_sp = solve_banded((1,1), ab, b) 
print("is our answer the same as scipy answer ({})".format(np.allclose(x, x_sp)))