Numpy '스마트'대칭 행렬

Numpy '스마트'대칭 행렬

---

쓸 [j][i]때 자동으로 (투명하게) 위치를 채우는 numpy에 스마트하고 공간 효율적인 대칭 행렬 [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

자동 Hermitian도 좋지만 글을 쓰는 시점에는 필요하지 않습니다.

---

계산을 수행하기 직전에 행렬을 대칭화할 여유가 있다면 다음이 상당히 빠릅니다.

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

이것은 합리적인 가정 하에서 작동합니다 (예 : 실행 전에 a[0, 1] = 42모순되는 것과 둘 다 수행하지 않음 ).a[1, 0] = 123symmetrize

투명 대칭 화가 정말로 필요한 경우 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…

---

numpy에서 대칭 행렬의 최적 처리에 대한보다 일반적인 문제도 저를 괴롭 혔습니다.

그것을 조사한 후 대답은 아마도 대칭 행렬에 대한 기본 BLAS 루틴이 지원하는 메모리 레이아웃에 의해 numpy가 다소 제한된다는 것입니다.

루틴이 할 몇 가지 BLAS는 대칭 행렬의 계산 속도를 대칭을 이용하지만, 그들은 여전히 전체 매트릭스, 같은 메모리 구조를 사용 n^2하기보다는 공간 n(n+1)/2. 그들은 행렬이 대칭이며 위쪽 또는 아래쪽 삼각형의 값만 사용하라는 말을 듣습니다.

일부 scipy.linalg루틴은 BLAS 루틴으로 전달되는 플래그 (예 : sym_pos=Trueon linalg.solve)를 허용 하지만 numpy에서 이에 대한 더 많은 지원이 좋을 것입니다. 특히 그램 행렬을 허용하는 DSYRK (대칭 순위 k 업데이트)와 같은 루틴에 대한 래퍼 dot (MT, M)보다 상당히 빠르게 계산됩니다.

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

---

대칭 행렬을 저장하는 잘 알려진 방법이 많이 있으므로 n ^ 2 저장 요소를 차지할 필요가 없습니다. 또한 이러한 수정 된 저장 수단에 액세스하기 위해 일반적인 작업을 다시 작성할 수 있습니다. 결정적인 작업은 Golub and 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은 저장된 V에 액세스하여 계산하는 방법을 보여주는 알고리즘 1.2.3을 제공합니다 y = V x + y.

Golub과 Van Loan은 또한 행렬을 대각선 우세한 형태로 저장하는 방법을 제공합니다. 이것은 스토리지를 저장하지 않지만 특정 다른 종류의 작업에 대한 준비된 액세스를 지원합니다.

---

이것은 평범한 파이썬이고 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

---

채워 [i][j]지면 파이썬 적으로 채우는 것은 사소한 일입니다 [j][i]. 저장 문제는 조금 더 흥미 롭습니다. 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)

```

참고 URL : https://stackoverflow.com/questions/2572916/numpy-smart-symmetric-matrix