Numpy 'smart'대칭 행렬

Question

[i][j]에 쓸 때 [j][i]에 자동으로 (그리고 투명하게) 위치를 채우는 스마트하고 공간 효율적인 대칭 행렬이 numpy에 있습니까?

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

, 다음이 합리적으로 빨리해야한다 :이 같은 symmetrize을 실행하기 전에 두 a[0, 1] = 42과 모순 a[1, 0] = 123을 수행하지 합리적인 가정 (아래 작동

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

).

당신이 정말 투명 대칭 화를해야하는 경우 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, :] 요소를 설정하는, a[:, 1] = -1 같은 더 복잡한 과제를 처리합니다. 코드가 약간 파이썬 3 실행하기 전에 적응해야한다, 그래서 파이썬 3, def …(…, (i, j),…)을 쓰기의 가능성을 제거하는 것이

참고 : def __setitem__(self, indexes, value): (i, j) = indexes ... NumPy와 대칭 행렬의 최적의 치료

Answer 2

더 많은 일반적인 문제는 나도 도청 .

그것으로보고 후에, 나는 대답은 아마 어느 정도 대칭 행렬에 대한 기본 BLAS 루틴 supportd 메모리 레이아웃에 의해 제한되는 NumPy와 생각합니다.

일부 BLAS 루틴은 대칭 행렬의 계산 속도를 높이기 위해 대칭을 사용하지만 전체 행렬과 동일한 메모리 구조, 즉 n(n+1)/2이 아닌 n^2 공간을 사용합니다. 그냥 그들은 행렬이 대칭 및 상하 한 삼각형 중 하나에서만 값을 사용하는 것을 말했다 얻을.

scipy.linalg 루틴의 일부를

는 플래그를 수용 할 NumPy와이에 대한 더 많은 지원이 좋은 것입니다 있지만 DSYRK (대칭 순위 K 업데이트)와 같은 루틴 특히 래퍼에, BLAS 루틴에 전달받을 ( linalg.solve에 sym_pos=True 등) 그램 매트릭스가 도트 (MT, M)보다 더 빨리 계산 될 수있게합니다.

(시간 및/또는 공간에 2x 상수 계수를 최적화하는 것에 대해 걱정할 필요는 없겠지만 단일 시스템에서 관리 할 수있는 문제의 크기에 대한 차이를 만들 수 있습니다 ...)

Answer 3

이 일반 파이썬과 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 and Van Loan, 매트릭스 계산, 3 판, 1996, Johns Hopkins University Press, 섹션 1.27-1.2.9. 예를 들어, 대칭 행렬에서 양식 (1.2.2)을 인용하면 i >= j에 A = [a_{i,j} ] 만 저장하면됩니다. 이어서, 상기 매트릭스를 들고 가정 벡터가 V로 표시되고있는 n 개의 별 N이다

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

이 1 인덱싱 가정에 a_{i,j} 넣어.

Golub 및 Van Loan은 저장된 V에 액세스하여 y = V x + y을 계산하는 방법을 보여주는 알고리즘 1.2.3을 제공합니다.

Golub 및 Van Loan은 또한 대각선 지배적 인 형태로 매트릭스를 저장하는 방법을 제공합니다. 이렇게하면 저장 용량을 절약 할 수는 있지만 특정 다른 종류의 작업에 대해서는 준비된 액세스를 지원합니다.

Answer 5

[j][i]이 채워지면 파이썬으로 [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) 

```접수로이 문제를 해결하는 경우