Самый быстрый способ создать разреженную матрицу вида A. T * diag(b) * A + C?
Я пытаюсь оптимизировать фрагмент кода, который решает большую разреженную нелинейную систему, используя метод внутренней точки. На этапе обновления это включает вычисление матрицы Гессена H, градиента g, а затем решение для d в H * d = -g, чтобы получить новое направление поиска.
Матрица Гессена имеет симметричную трехдиагональную структуру вида:
A. T * diag (b) * A + C
Я убежал.line_profiler о конкретной функции в вопросе:
Line # Hits Time Per Hit % Time Line Contents
==================================================
386 def _direction(n, res, M, Hsig, scale_var, grad_lnprior, z, fac):
387
388 # gradient
389 44 1241715 28220.8 3.7 g = 2 * scale_var * res - grad_lnprior + z * np.dot(M.T, 1. / n)
390
391 # hessian
392 44 3103117 70525.4 9.3 N = sparse.diags(1. / n ** 2, 0, format=FMT, dtype=DTYPE)
393 44 18814307 427597.9 56.2 H = - Hsig - z * np.dot(M.T, np.dot(N, M)) # slow!
394
395 # update direction
396 44 10329556 234762.6 30.8 d, fac = my_solver(H, -g, fac)
397
398 44 111 2.5 0.0 return d, fac
Глядя на результат, становится ясно, что построение
H является, безусловно, самым дорогостоящим шагом - оно занимает значительно больше времени, чем собственно решение для нового направления. Hsig и M являются разреженными матрицами CSC, n - плотный вектор и z - скаляр. Решатель, который я использую, требует, чтобы H был либо CSC, либо CSR разреженной матрицей.
Вот функция, которая производит некоторые игрушечные данные с теми же форматами, размерами и разреженностью, что и мои вещественные матрицы:
import numpy as np
from scipy import sparse
def make_toy_data(nt=200000, nc=10):
d0 = np.random.randn(nc * (nt - 1))
d1 = np.random.randn(nc * (nt - 1))
M = sparse.diags((d0, d1), (0, nc), shape=(nc * (nt - 1), nc * nt),
format='csc', dtype=np.float64)
d0 = np.random.randn(nc * nt)
Hsig = sparse.diags(d0, 0, shape=(nc * nt, nc * nt), format='csc',
dtype=np.float64)
n = np.random.randn(nc * (nt - 1))
z = np.random.randn()
return Hsig, M, n, z
И вот мой оригинальный подход к построению H:
def original(Hsig, M, n, z):
N = sparse.diags(1. / n ** 2, 0, format='csc')
H = - Hsig - z * np.dot(M.T, np.dot(N, M)) # slow!
return H
Время:
%timeit original(Hsig, M, n, z)
# 1 loops, best of 3: 483 ms per loop
Есть ли более быстрый способ построить эту матрицу?
2 ответов:
Я приближаюсь к 4-кратному ускорению вычисления произведения
M.T * D * Mиз трех диагональных массивов. Еслиd0иd1являются главной и верхней диагональюM, аdявляется главной диагональюD, то следующий код создаетM.T * D * Mнепосредственно:def make_tridi_bis(d0, d1, d, nc=10): d00 = d0*d0*d d11 = d1*d1*d d01 = d0*d1*d len_ = d0.size data = np.empty((3*len_ + nc,)) indices = np.empty((3*len_ + nc,), dtype=np.int) # Fill main diagonal data[:2*nc:2] = d00[:nc] indices[:2*nc:2] = np.arange(nc) data[2*nc+1:-2*nc:3] = d00[nc:] + d11[:-nc] indices[2*nc+1:-2*nc:3] = np.arange(nc, len_) data[-2*nc+1::2] = d11[-nc:] indices[-2*nc+1::2] = np.arange(len_, len_ + nc) # Fill top diagonal data[1:2*nc:2] = d01[:nc] indices[1:2*nc:2] = np.arange(nc, 2*nc) data[2*nc+2:-2*nc:3] = d01[nc:] indices[2*nc+2:-2*nc:3] = np.arange(2*nc, len_+nc) # Fill bottom diagonal data[2*nc:-2*nc:3] = d01[:-nc] indices[2*nc:-2*nc:3] = np.arange(len_ - nc) data[-2*nc::2] = d01[-nc:] indices[-2*nc::2] = np.arange(len_ - nc ,len_) indptr = np.empty((len_ + nc + 1,), dtype=np.int) indptr[0] = 0 indptr[1:nc+1] = 2 indptr[nc+1:len_+1] = 3 indptr[-nc:] = 2 np.cumsum(indptr, out=indptr) return sparse.csr_matrix((data, indices, indptr), shape=(len_+nc, len_+nc))Если ваша матрица
Mбыла в формате CSR, вы можете извлечьd0иd1какd0 = M.data[::2]иd1 = M.data[1::2], я модифицировал вашу игрушечную процедуру создания данных, чтобы вернуть эти массивы, и вот что я получаю:In [90]: np.allclose((M.T * sparse.diags(d, 0) * M).A, make_tridi_bis(d0, d1, d).A) Out[90]: True In [92]: %timeit make_tridi_bis(d0, d1, d) 10 loops, best of 3: 124 ms per loop In [93]: %timeit M.T * sparse.diags(d, 0) * M 1 loops, best of 3: 501 ms per loop
Целое цель приведенного выше кода состоит в том, чтобы воспользоваться структурой ненулевых записей. Если вы нарисуете диаграмму матриц, которые вы умножаете вместе, относительно легко убедить себя, что главные (
d_0) и верхние и нижние (d_1) диагонали результирующей трехдиагональной матрицы просто:Остальная часть кода в этой функции просто строит непосредственно трехдиагональную матрицу, так как вызовd_0 = np.zeros((len_ + nc,)) d_0[:len_] = d00 d_0[-len_:] += d11 d_1 = d01sparse.diagsс указанными выше данными в несколько раз медленнее.
Я попробовал запустить ваш тестовый случай и у меня возникли проблемы с
np.dot(N, M). Я не копался в этом, но я думаю, что у моего numpy / sparse combo (оба довольно новые) были проблемы с использованиемnp.dotна разреженных массивах.Но
H = -Hsig - z*M.T.dot(N.dot(M))работает просто отлично. При этом используетсяsparse dot.Я не запускал профиль, но вот тайминги Ipython для нескольких частей. Для генерации данных требуется больше времени, чем для этой двойной точки.
In [37]: timeit Hsig,M,n,z=make_toy_data() 1 loops, best of 3: 2 s per loop In [38]: timeit N = sparse.diags(1. / n ** 2, 0, format='csc') 1 loops, best of 3: 377 ms per loop In [39]: timeit H = -Hsig - z*M.T.dot(N.dot(M)) 1 loops, best of 3: 1.55 s per loop
Hявляется<2000000x2000000 sparse matrix of type '<type 'numpy.float64'>' with 5999980 stored elements in Compressed Sparse Column format>
Comments