Numpy'smart'symmetric matrix

Question

numpy中是否有一個智能且節省空間的對稱矩陣，當寫入[i][j]時，會自動（並透明地）填充[j][i]處的位置？

import numpy 
a = numpy.symmetric((3, 3)) 
a[0][1] = 1 
a[1][0] == a[0][1] 
# True 
print(a) 
# [[0 1 0], [1 0 0], [0 0 0]] 

assert numpy.all(a == a.T) # for any symmetric matrix 

一個自動Hermitian也會很好，雖然我不會在寫作時需要它。

Answer 1

如果你能負擔得起只是在做計算之前，對稱化矩陣，下面應該是相當快的：

def symmetrize(a): 
    return a + a.T - numpy.diag(a.diagonal()) 

這個工程在合理的假設（如運行symmetrize之前沒有做既a[0, 1] = 42和矛盾a[1, 0] = 123 ）。

如果你真的需要一個透明對稱，你可能會考慮繼承numpy.ndarray，只是重新定義__setitem__：

class SymNDArray(numpy.ndarray): 
    def __setitem__(self, (i, j), value): 
     super(SymNDArray, self).__setitem__((i, j), value)      
     super(SymNDArray, self).__setitem__((j, i), value)      

def symarray(input_array): 
    """ 
    Returns a symmetrized version of the array-like input_array. 
    Further assignments to the array are automatically symmetrized. 
    """ 
    return symmetrize(numpy.asarray(input_array)).view(SymNDArray) 

# Example: 
a = symarray(numpy.zeros((3, 3))) 
a[0, 1] = 42 
print a # a[1, 0] == 42 too! 

（或矩陣來代替陣列等效，根據您的需要）。這種方法甚至可以處理更復雜的分配，如a[:, 1] = -1，它可以正確設置a[1, :]元素。

注意，Python 3中取出寫def …(…, (i, j),…)的可能性，所以代碼必須與Python 3運行前會略有調整：def __setitem__(self, indexes, value): (i, j) = indexes ...

Answer 2

在numpy的對稱矩陣的最佳治療方案的更普遍的問題竊聽我太。

尋找到它之後，我想答案是可能是numpy的是有些由底層BLAS例程對稱矩陣supportd內存佈局的限制。

儘管一些BLAS例程利用對稱性來加速對稱矩陣上的計算，但它們仍然使用與全矩陣相同的存儲器結構，即n^2空間而不是n(n+1)/2。只是他們被告知矩陣是對稱的，並且只使用上三角或下三角中的值。

一些scipy.linalg例程不接受標誌（像sym_pos=True上linalg.solve），其中獲得傳遞給BLAS例程，雖然這在numpy的更多的支持將是很好，尤其是包裝的程序像DSYRK（對稱k級更新） ，這將允許計算一個格點矩陣比點（MT，M）快得多。

（可能看起來很挑剔，擔心在時間和/或空間上優化2倍恆定因子，但它可以在單個機器上管理多大問題的閾值上有所不同......）

Answer 3

這是普通的Python，而不是numpy的，但我只是扔在一起例行填補 對稱矩陣（和測試程序，以確保它是正確的）：

import random 

# fill a symmetric matrix with costs (i.e. m[x][y] == m[y][x] 
# For demonstration purposes, this routine connect each node to all the others 
# Since a matrix stores the costs, numbers are used to represent the nodes 
# so the row and column indices can represent nodes 

def fillCostMatrix(dim):  # square array of arrays 
    # Create zero matrix 
    new_square = [[0 for row in range(dim)] for col in range(dim)] 
    # fill in main diagonal 
    for v in range(0,dim): 
     new_square[v][v] = random.randrange(1,10) 

    # fill upper and lower triangles symmetrically by replicating diagonally 
    for v in range(1,dim): 
     iterations = dim - v 
     x = v 
     y = 0 
     while iterations > 0: 
      new_square[x][y] = new_square[y][x] = random.randrange(1,10) 
      x += 1 
      y += 1 
      iterations -= 1 
    return new_square 

# sanity test 
def test_symmetry(square): 
    dim = len(square[0]) 
    isSymmetric = '' 
    for x in range(0, dim): 
     for y in range(0, dim): 
      if square[x][y] != square[y][x]: 
       isSymmetric = 'NOT' 
    print "Matrix is", isSymmetric, "symmetric" 

def showSquare(square): 
    # Print out square matrix 
    columnHeader = ' ' 
    for i in range(len(square)): 
     columnHeader += ' ' + str(i) 
    print columnHeader 

    i = 0; 
    for col in square: 
     print i, col # print row number and data 
     i += 1 

def myMain(argv): 
    if len(argv) == 1: 
     nodeCount = 6 
    else: 
     try: 
      nodeCount = int(argv[1]) 
     except: 
      print "argument must be numeric" 
      quit() 

    # keep nodeCount <= 9 to keep the cost matrix pretty 
    costMatrix = fillCostMatrix(nodeCount) 
    print "Cost Matrix" 
    showSquare(costMatrix) 
    test_symmetry(costMatrix) # sanity test 
if __name__ == "__main__": 
    import sys 
    myMain(sys.argv) 

# vim:tabstop=8:shiftwidth=4:expandtab 

Answer 4

有許多的好已知的存儲對稱矩陣的方法使得它們不需要佔用n^2個存儲元件。此外，重寫通用操作來訪問這些修改的存儲方法是可行的。最終的工作是Golub和Van Loan，Matrix Computations，1996年第3版，Johns Hopkins University Press，第1.27-1.2.9節。例如，從形式（1.2.2）引用它們，在對稱矩陣中，只需要存儲A = [a_{i,j} ]就可以得到i >= j。然後，假設矢量保持矩陣被表示爲V，並且A爲n乘n，把a_{i,j}在

V[(j-1)n - j(j-1)/2 + i] 

這假定1索引。

Golub和Van Loan提供了一個算法1.2.3，它顯示瞭如何訪問這樣存儲的V來計算y = V x + y。

Golub和Van Loan也提供了一種以對角佔優形式存儲矩陣的方法。這不會節省存儲空間，但支持對某些其他類型的操作進行訪問。

Answer 5

如果填寫了[j][i]，那麼通過Python填充[i][j]是微不足道的。存儲問題更有趣一點。可以使用packed屬性來擴充numpy數組類，這對於保存存儲和稍後讀取數據都很有用。

class Sym(np.ndarray): 

    # wrapper class for numpy array for symmetric matrices. New attribute can pack matrix to optimize storage. 
    # Usage: 
    # If you have a symmetric matrix A as a shape (n,n) numpy ndarray, Sym(A).packed is a shape (n(n+1)/2,) numpy array 
    # that is a packed version of A. To convert it back, just wrap the flat list in Sym(). Note that Sym(Sym(A).packed) 


    def __new__(cls, input_array): 
     obj = np.asarray(input_array).view(cls) 

     if len(obj.shape) == 1: 
      l = obj.copy() 
      p = obj.copy() 
      m = int((np.sqrt(8 * len(obj) + 1) - 1)/2) 
      sqrt_m = np.sqrt(m) 

      if np.isclose(sqrt_m, np.round(sqrt_m)): 
       A = np.zeros((m, m)) 
       for i in range(m): 
        A[i, i:] = l[:(m-i)] 
        A[i:, i] = l[:(m-i)] 
        l = l[(m-i):] 
       obj = np.asarray(A).view(cls) 
       obj.packed = p 

      else: 
       raise ValueError('One dimensional input length must be a triangular number.') 

     elif len(obj.shape) == 2: 
      if obj.shape[0] != obj.shape[1]: 
       raise ValueError('Two dimensional input must be a square matrix.') 
      packed_out = [] 
      for i in range(obj.shape[0]): 
       packed_out.append(obj[i, i:]) 
      obj.packed = np.concatenate(packed_out) 

     else: 
      raise ValueError('Input array must be 1 or 2 dimensional.') 

     return obj 

    def __array_finalize__(self, obj): 
     if obj is None: return 
     self.packed = getattr(obj, 'packed', None) 

```