Умная симметричная матрица Numpy

Ответ 1

Если вы можете позволить себе симметризировать матрицу непосредственно перед выполнением вычислений, следующее должно быть достаточно быстрым:

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

Это работает при разумных предположениях (например, не делая как a[0, 1] = 42, так и противоречивый a[1, 0] = 123 перед запуском symmetrize).

Если вам действительно нужна прозрачная симметризация, вы можете рассмотреть подклассификацию 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 код

Ответ 2

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

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

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

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

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

Ответ 3

Существует ряд хорошо известных способов хранения симметричных матриц, поэтому им не нужно занимать n ^ 2 элемента хранения. Кроме того, можно переписать общие операции для доступа к этим пересмотренным средствам хранения. Окончательной работой являются Голуб и Ван Лоан, Matrix Computations, 3-е издание 1996 года, Университетский университет Джона Хопкинса, разделы 1.27-1.2.9. Например, цитируя их из формы (1.2.2), в симметричной матрице нужно хранить A = [a_{i,j} ] для i >= j. Тогда, предполагая, что вектор, удерживающий матрицу, обозначается V, а A - n-by-n, поместите a_{i,j} в

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

Это предполагает 1-индексирование.

Голуб и Ван Лоан предлагают Алгоритм 1.2.3, который показывает, как получить доступ к такой сохраненной V, чтобы вычислить y = V x + y.

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

Ответ 4

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

Ответ 5

Это тривиально для 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)

`` `