Включает в себя "умные" симметричной матрицы

более общая проблема оптимальной обработки симметричных матриц в numpy тоже меня беспокоила.

после изучения этого я думаю, что ответ, вероятно, заключается в том, что numpy несколько ограничен поддержкой макета памяти базовыми подпрограммами BLAS для симметричных матриц.

хотя некоторые процедуры BLAS используют симметрию для ускорения вычислений на симметричных матрицах, они по-прежнему используют ту же структуру памяти, что и полная матрица, то есть n^2 пространства, а не n(n+1)/2. Просто им говорят, что матрица симметрична и использовать только значения в верхнем или нижнем треугольнике.

часть scipy.linalg подпрограммы принимают флаги (например,sym_pos=True on linalg.solve), которые переходят к подпрограммам BLAS, хотя больше поддержки для этого в numpy было бы неплохо, в частности, обертки для таких подпрограмм, как DSYRK (обновление симметричного ранга k), что позволило бы вычислять матрицу Грама намного быстрее, чем dot(M. T, M).

(кажется nitpicky беспокоиться об оптимизации для постоянного фактора 2x по времени и / или пространству, но это может иметь значение для этого порога того, насколько большой проблемой вы можете управлять на одной машине...)

Существует ряд известных способов хранения симметричных матриц, поэтому им не нужно занимать N^2 элементов хранения. Кроме того, можно переписать общие операции для доступа к этим пересмотренным средствам хранения. Окончательная работа-Голуб и Ван займ,Матричные Вычисления, 3-е издание 1996 года, издательство Университета Джонса Хопкинса, разделы 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.

Голуб и Ван займ также обеспечивают способ хранения матрицы в диагональной доминирующей форме. Это не экономит память, но поддерживает готовый доступ для некоторых других видов операций.

Это простой 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

это тривиально, чтобы Pythonically заполнить [i][j] если [j][i] заполняется. Вопрос хранения немного интереснее. Можно увеличить класс массива numpy с помощью packed атрибут, который полезен как для сохранения памяти, так и для последующего чтения данных.

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)

"'