Обычно я бы не стал бить дохлую лошадь, но бывает, что мой не векторизованный подход к решению вашего другого вопроса имеет некоторые преимущества, когда дела становятся большими. Поскольку я фактически заполнял матрицу коэффициентов по одному элементу за раз, ее очень легко перевести в формат разреженной матрицы COO, который можно эффективно преобразовать в CSC и решить. Это делает следующее:
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)
Это, конечно, гораздо более многословно, чем построение из векторизованных блоков, как это делает Theodoros Zelleke, но происходит интересная вещь, когда вы учитываете оба подхода:
![введите здесь описание изображения](https://i.stack.imgur.com/QJ5CY.png)
Во-первых, и это (очень) приятно, время линейно масштабируется в обоих решениях, как и следовало ожидать при использовании разреженного подхода. Но решение, которое я дал в этом ответе, всегда быстрее, особенно для больших n
. Просто для удовольствия я также рассчитал плотный подход Теодорос Зеллеке из другого вопроса, который дает этот хороший график, показывающий различное масштабирование обоих типов решений и то, как очень рано, где-то около n = 75
, решение здесь должно быть вашим выбором:
![введите здесь описание изображения](https://i.stack.imgur.com/8HAxc.png)
Я недостаточно знаю о scipy.sparse
, чтобы действительно понять, почему существуют различия между двумя разреженными подходами, хотя я сильно подозреваю об использовании разреженных матриц формата LIL. Может быть очень незначительный прирост производительности, хотя в коде много компактности, если перевести ответ TheodorosZelleke в формат COO. Но это остается в качестве упражнения для ОП!
person
Jaime
schedule
17.01.2013