Ma trận đối xứng 'thông minh' Numpy

Question

Có ma trận đối xứng thông minh và không gian hiệu quả trong khối ô vuông tự động (và trong suốt) lấp đầy vị trí tại [j][i] khi viết [i][j]?

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 

Một Hermitian tự động cũng sẽ tốt, mặc dù tôi không cần nó tại thời điểm viết.

Answer 1

Nếu bạn có thể đủ khả năng để làm cho cân đối các ma trận ngay trước khi làm phép tính, sau đây nên được hợp lý nhanh:

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

này hoạt động theo các giả định hợp lý (chẳng hạn như không làm cả hai a[0, 1] = 42 và mâu thuẫn a[1, 0] = 123 trước khi chạy symmetrize).

Nếu bạn thực sự cần một symmetrization trong suốt, bạn có thể xem xét subclassing numpy.ndarray và chỉ cần xác định lại __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! 

(hoặc tương đương với ma trận thay vì mảng, tùy thuộc vào nhu cầu của bạn). Cách tiếp cận này thậm chí còn xử lý các bài tập phức tạp hơn, chẳng hạn như a[:, 1] = -1, cách đặt chính xác các yếu tố a[1, :].

Lưu ý rằng Python 3 loại bỏ khả năng viết def …(…, (i, j),…), vì vậy các mã có được một chút thích nghi trước khi chạy với Python 3: def __setitem__(self, indexes, value): (i, j) = indexes ...

Answer 2

Các vấn đề chung của điều trị tối ưu của ma trận đối xứng trong NumPy nghe trộm tôi quá .

Sau khi nhìn vào nó, tôi nghĩ câu trả lời có lẽ là một phần nào đó bị hạn chế bởi sự hỗ trợ bố trí bộ nhớ bởi các thói quen BLAS cơ bản cho các ma trận đối xứng.

Trong khi một số quy trình BLAS khai thác đối xứng để tăng tốc độ tính toán trên ma trận đối xứng, chúng vẫn sử dụng cùng cấu trúc bộ nhớ như ma trận đầy đủ, nghĩa là, n^2 không gian thay vì n(n+1)/2. Chỉ cần họ được cho biết rằng ma trận là đối xứng và chỉ sử dụng các giá trị trong tam giác trên hoặc dưới.

Một số scipy.linalg thói quen làm chấp nhận lá cờ (như sym_pos=True trên linalg.solve) mà có được thông qua vào thói quen BLAS, mặc dù nhiều hỗ trợ cho điều này trong NumPy sẽ được tốt đẹp, trong giấy gói đặc biệt đối với thói quen như DSYRK (cấp bậc đối xứng k cập nhật) , cho phép ma trận Gram được tính toán một chút công bằng nhanh hơn dấu chấm (MT, M).

(Có vẻ như nitpicky phải lo lắng về việc tối ưu hóa hệ số không đổi 2x về thời gian và/hoặc không gian, nhưng nó có thể tạo ra sự khác biệt cho mức độ lớn mà bạn có thể quản lý trên một máy ...)

Answer 3

Đây là python đơn giản và không numPy, nhưng tôi chỉ ném cùng một thói quen để điền một ma trận đối xứng (và một chương trình thử nghiệm để chắc chắn rằng nó là chính xác):

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

có một số tốt các cách lưu trữ ma trận đối xứng để chúng không cần chiếm giữ^^ yếu tố lưu trữ. Hơn nữa, nó là khả thi để viết lại các hoạt động phổ biến để truy cập các phương tiện lưu trữ sửa đổi.Tác phẩm cuối cùng là Golub và Văn Loan, Tính toán ma trận, ấn bản lần thứ 3 năm 1996, Nhà xuất bản Đại học Johns Hopkins, phần 1.27-1.2.9. Ví dụ, trích dẫn chúng từ biểu mẫu (1.2.2), trong ma trận đối xứng chỉ cần lưu trữ A = [a_{i,j} ] cho i >= j. Sau đó, giả sử vector giữ ma trận được ký hiệu V, và rằng A là n-by-n, đặt a_{i,j} trong

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

này giả định 1-lập chỉ mục.

Golub và Van Loan cung cấp Thuật toán 1.2.3 cho biết cách truy cập V được lưu trữ như vậy để tính y = V x + y.

Golub và Van Loan cũng cung cấp cách lưu trữ ma trận ở dạng chi phối chéo. Điều này không lưu trữ, nhưng hỗ trợ truy cập sẵn sàng cho một số loại hoạt động khác.

Answer 5

Điều nhỏ nhặt là điền vào [i][j] nếu điền [j][i]. Câu hỏi lưu trữ thú vị hơn một chút. Người ta có thể tăng thêm các lớp mảng numpy với một thuộc tính packed hữu ích cả để lưu trữ và sau đó đọc dữ liệu.

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) 

`` `