Używanie scipy rzadkich macierzy do rozwiązania układu równań

Question

Jest to kontynuacja do How to set up and solve simultaneous equations in python, ale mam wrażenie, że zasługuje na własną reputację za jakąkolwiek odpowiedź.

Dla stałej liczby całkowitej n, mam zestaw równoczesnych równań 2(n-1) w następujący sposób.

M(p) = 1+((n-p-1)/n)*M(n-1) + (2/n)*N(p-1) + ((p-1)/n)*M(p-1) 

N(p) = 1+((n-p-1)/n)*M(n-1) + (p/n)*N(p-1) 

M(1) = 1+((n-2)/n)*M(n-1) + (2/n)*N(0) 

N(0) = 1+((n-1)/n)*M(n-1) 

M(p) jest zdefiniowany 1 <= p <= n-1. N(p) jest zdefiniowany dla 0 <= p <= n-2. Zauważ również, że p jest po prostu stałą liczbą całkowitą w każdym równaniu, więc cały system jest liniowy.

Podano kilka bardzo miłych odpowiedzi dotyczących konfiguracji układu równań w pytonie. Jednak system jest rzadki i chciałbym go rozwiązać dla dużych n. Jak mogę zamiast tego użyć reprezentacji rzadkiej macierzy scipy i na przykład http://docs.scipy.org/doc/scipy/reference/sparse.linalg.html?

Answer 1

nie będę normalnie zachować pokonując martwego konia, ale zdarza się, że mój non-wektorowy podejście do rozwiązywania inne pytanie, ma pewne zasługi, gdy sprawy być wielkim. Ponieważ faktycznie wypełniałem macierz współczynników po jednym elemencie na raz, bardzo łatwo jest przetłumaczyć na rzadki format macierzy COO, który można skutecznie przekształcić w CSC i rozwiązać. Dodaje się robi:

import scipy.sparse 

def sps_solve(n) : 
    # Solution vector is [N[0], N[1], ..., N[n - 2], M[1], M[2], ..., M[n - 1]] 
    n_pos = lambda p : p 
    m_pos = lambda p : p + n - 2 
    data = [] 
    row = [] 
    col = [] 
    # p = 0 
    # n * N[0] + (1 - n) * M[n-1] = n 
    row += [n_pos(0), n_pos(0)] 
    col += [n_pos(0), m_pos(n - 1)] 
    data += [n, 1 - n] 
    for p in xrange(1, n - 1) : 
     # n * M[p] + (1 + p - n) * M[n - 1] - 2 * N[p - 1] + 
     # (1 - p) * M[p - 1] = n 
     row += [m_pos(p)] * (4 if p > 1 else 3) 
     col += ([m_pos(p), m_pos(n - 1), n_pos(p - 1)] + 
       ([m_pos(p - 1)] if p > 1 else [])) 
     data += [n, 1 + p - n , -2] + ([1 - p] if p > 1 else []) 
     # n * N[p] + (1 + p -n) * M[n - 1] - p * N[p - 1] = n 
     row += [n_pos(p)] * 3 
     col += [n_pos(p), m_pos(n - 1), n_pos(p - 1)] 
     data += [n, 1 + p - n, -p] 
    if n > 2 : 
     # p = n - 1 
     # n * M[n - 1] - 2 * N[n - 2] + (2 - n) * M[n - 2] = n 
     row += [m_pos(n-1)] * 3 
     col += [m_pos(n - 1), n_pos(n - 2), m_pos(n - 2)] 
     data += [n, -2, 2 - n] 
    else : 
     # p = 1 
     # n * M[1] - 2 * N[0] = n 
     row += [m_pos(n - 1)] * 2 
     col += [m_pos(n - 1), n_pos(n - 2)] 
     data += [n, -2] 
    coeff_mat = scipy.sparse.coo_matrix((data, (row, col))).tocsc() 
    return scipy.sparse.linalg.spsolve(coeff_mat, 
             np.ones(2 * (n - 1)) * n) 

Jest oczywiście dużo bardziej gadatliwy niż budowanie go od vectorized bloków, tak jak TheodorosZelleke, ale ciekawa rzecz dzieje się, gdy czas obu podejść:

Po pierwsze, a to jest (bardzo) miłe, czas skaluje się liniowo w obu rozwiązaniach, tak jak można by się spodziewać, stosując rzadkie podejście. Ale rozwiązanie, które dałem w tej odpowiedzi, jest zawsze szybsze, bardziej dla większych n s. Dla samej zabawy, odmieniłem także gęste podejście TheodorosZelleke'a z drugiego pytania, które daje ten przyjemny wykres pokazujący różne skalowanie obu typów rozwiązań i jak bardzo wcześnie, gdzieś w okolicach n = 75, rozwiązanie tutaj powinno być twoim wyborem:

nie wiem wystarczająco dużo o scipy.sparse aby naprawdę zrozumieć, dlaczego różnice między tymi dwoma podejściami rzadkie, chociaż podejrzewam, że ciężko z wykorzystaniem formatu LIL rozrzedzony matryc. Może pojawić się niewielki wzrost wydajności, pomimo dużej zwartości kodu, poprzez przekształcenie odpowiedzi TheodorosZelleke na format COO. Ale to jest pozostawione jako ćwiczenie dla OP!

Answer 2

Jest pewien kod, na który patrzyłem wcześniej: http://jkwiens.com/heat-equation-using-finite-difference/ Jego funkcja implementuje metodę różnic skończonych w celu rozwiązania równania ciepła przy użyciu pakietu scipy sparse matrix.

Answer 3

To jest rozwiązanie przy użyciu scipy.sparse. Niestety problem nie został tutaj opisany. Aby zrozumieć to rozwiązanie, przyszli użytkownicy muszą najpierw sprawdzić problem pod linkiem podanym w pytaniu.

rozwiązanie wykorzystujące scipy.sparse:

from scipy.sparse import spdiags, lil_matrix, vstack, hstack 
from scipy.sparse.linalg import spsolve 
import numpy as np 


def solve(n): 
    nrange = np.arange(n) 
    diag = np.ones(n-1) 

    # upper left block 
    n_to_M = spdiags(-2. * diag, 0, n-1, n-1) 

    # lower left block 
    n_to_N = spdiags([n * diag, -nrange[-1:0:-1]], [0, 1], n-1, n-1) 

    # upper right block 
    m_to_M = lil_matrix(n_to_N) 
    m_to_M[1:, 0] = -nrange[1:-1].reshape((n-2, 1)) 

    # lower right block 
    m_to_N = lil_matrix((n-1, n-1)) 
    m_to_N[:, 0] = -nrange[1:].reshape((n-1, 1)) 

    # build A, combine all blocks 
    coeff_mat = hstack(
         (vstack((n_to_M, n_to_N)), 
         vstack((m_to_M, m_to_N)))) 

    # const vector, right side of eq. 
    const = n * np.ones((2 * (n-1),1)) 

    return spsolve(coeff_mat.tocsr(), const).reshape((-1,1))