Numpy 'smart' symetryczna matryca

Question

Czy istnieje inteligentna i efektywna przestrzennie symetryczna macierz w numpy, która automatycznie (i przeźroczyście) wypełnia pozycję na [j][i], gdy zostanie zapisana [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 

Również automatyczny heritianizm byłby niezły, chociaż nie będę tego potrzebował w chwili pisania.

Answer 1

Jeśli możesz sobie pozwolić na symetryzować matrycę tylko przed wykonaniem obliczeń, następujące dokumenty powinny być w miarę szybko:

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

Działa to na rozsądnych założeniach (jak nie robi zarówno a[0, 1] = 42 i sprzeczne a[1, 0] = 123 przed uruchomieniem symmetrize).

Jeśli naprawdę potrzebują przejrzystych symetryczne, można rozważyć instacji numpy.ndarray i po prostu przedefiniowanie __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! 

(lub równowartość w matrycach zamiast tablic, w zależności od potrzeb). Takie podejście obsługuje nawet bardziej skomplikowane zadania, takie jak a[:, 1] = -1, które poprawnie ustawia elementy a[1, :].

Zauważ, że Python 3 usunięto możliwość pisania def …(…, (i, j),…), więc kod musi być nieznacznie dostosowane przed uruchomieniem z Python 3: def __setitem__(self, indexes, value): (i, j) = indexes ...

Answer 2

Im bardziej ogólny problem optymalnego leczenia macierzy symetrycznych w numpy podsłuch mnie zbyt .

Po przeanalizowaniu tego, myślę, że odpowiedź jest prawdopodobnie taka, że ​​numpy jest nieco ograniczony przez układ pamięci obsługiwany przez podstawowe procedury BLAS dla macierzy symetrycznych.

Podczas gdy niektóre procedury BLAS wykorzystują symetrię do przyspieszenia obliczeń na macierzach symetrycznych, nadal używają tej samej struktury pamięci, co pełna macierz, czyli przestrzeń n^2, a nie n(n+1)/2. Dowiadują się, że macierz jest symetryczna i używa tylko wartości w górnym lub dolnym trójkącie.

Niektóre scipy.linalg rutyny akceptuje flagi (jak sym_pos=True na linalg.solve), który przejdzie do procedury Blas, chociaż większe wsparcie dla tego w numpy byłoby miło, w szczególności opakowań do rutyny jak DSYRK (Rank symetryczny aktualizacji k) , które pozwoliłyby na obliczenie macierzy Grama szybciej niż kropka (MT, M).

(Może wydawać się dziwaczny, aby martwić się o optymalizację w celu uzyskania 2-krotnego stałego współczynnika na czas i/lub przestrzeń, ale może mieć wpływ na to, jak duży problem można sobie poradzić na jednej maszynie ...)

Answer 3

jest to zwykły python i nie numpy, ale ja po prostu wyrzucił razem rutynowych wypełnić macierz symetryczną (i program testowy, aby upewnić się, że jest poprawne):

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

Istnieje wiele dobrze - znane sposoby przechowywania macierzy symetrycznych, dzięki czemu nie muszą zajmować n^2 elementów pamięci. Ponadto możliwe jest przepisanie typowych operacji w celu uzyskania dostępu do zmienionych sposobów przechowywania.Ostateczna praca to Golub and Van Loan, Matrix Computations, 3rd edition 1996, Johns Hopkins University Press, sekcje 1.27-1.2.9. Na przykład, cytując je z formularza (1.2.2), w macierzy symetrycznej wystarczy przechowywać A = [a_{i,j} ] dla i >= j. Następnie, zakładając, że wektorposiadających matrycę oznaczamy V, a A oznacza N, przez N, umieścić a_{i,j} w

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

Zakłada 1 indeksowanie.

Golub i Van Loan oferują algorytm 1.2.3, który pokazuje, jak uzyskać dostęp do takiego zapisanego V, aby obliczyć y = V x + y.

Golub i Van Loan zapewniają również sposób przechowywania matrycy w dominującej formie przekątnej. Nie oszczędza to pamięci, ale obsługuje łatwy dostęp do niektórych innych operacji.

Answer 5

Pythonically wypełnićjest trywialnie wypełnione [i][j] jeśli wypełni się [j][i]. Pytanie o przechowywanie jest trochę bardziej interesujące. Można wzmocnić numpy tablicę klasy atrybutem packed, która jest przydatna zarówno do zapisywania pamięci, jak i do późniejszego odczytu danych.

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) 

`` `