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

Question

Существует ли интеллектуальная и компактная симметричная матрица в numpy, которая автоматически (и прозрачно) заполняет позицию в [j][i] при записи в [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

Автоматический эрмитов тоже подойдет, хотя мне это и не понадобится во время написания.

Answer 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: def __setitem__(self, indexes, value): (i, j) = indexes…

Answer 2

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

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

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

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

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

Answer 3

Существует ряд известных способов хранения симметричных матриц, поэтому им не нужно занимать 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-индексирование.

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

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

Answer 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

Answer 5

Заполнение [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)

`` `