Мне нужно инвертировать большую плотную матрицу, которую я надеялся использовать для работы gmres Scipy. К счастью, плотная матрица A следует шаблону, и мне не нужно хранить матрицу в памяти. Класс LinearOperator позволяет нам построить объект, который действует как матрица для GMRES и может вычислить непосредственно векторное векторное произведение A*v. То есть, мы пишем функцию mv(v), которая берет на вход вектор v и возвращает mv(v) = A*v. Затем мы можем использовать класс LinearOperator для создания A_LinOp = LinearOperator(shape = shape, matvec = mv). Мы можем поместить линейный оператор в команду Scipy gmres для оценки матричных векторных продуктов без необходимости полностью загружать A в память.Scipy LinearOperator с несколькими входами
Документация для LinearOperator находится здесь: LinearOperator Documentation.
Вот моя проблема: написать подпрограмму для вычисления векторного векторного произведения mv(v) = A*v, мне нужен еще один вектор ввода C. Записи в A имеют форму A[i,j] = f(C[i] - C[j]). Итак, я действительно хочу, чтобы mv состоял из двух входов, один фиксированный векторный вход C и один вход переменной v, который мы хотим вычислить A*v.
MATLAB имеет подобную установку, где бы написать x = gmres(@(v) mv(v,C),b) где b правая рука сторона проблемы Ax = b, и mv это функция, которая принимает в качестве переменной ввода v, который мы хотим вычислить A*v и C является фиксированной, который нам нужен для сборки A.
Моя проблема заключается в том, что я не могу понять, как разрешить классу LinearOperator принимать два входа, одну переменную и одну «фиксированную», как я могу в MATLAB.
Есть ли способ сделать аналогичную операцию в SciPy? В качестве альтернативы, если кто-нибудь знает лучший способ обращения большой, плотной матрицы (50000, 50000), где записи следуют шаблону, я был бы очень признателен за любые предложения.
Спасибо!
EDIT: Я должен был сообщить эту информацию на самом деле. Матрица на самом деле (в виде блока) [A C; C^T 0], где находится AN x N (N большой) и C является N x 3, а 0 является 3 x 3 и C^T является транспонированной C. Этот массив C представляет собой тот же массив, что и упомянутый выше. Записи A следуют рисунку A[i,j] = f(C[i] - C[j]).
Я написал mv(v,C) идти по строкам построить A*v[i] для i=0,N, вычисляя сумму f(C[i]-C[j)*v[j] (на самом деле, я numpy.dot(FC,v) где FC[j] = f(C[i]-C[j]), который хорошо работает). Затем в конце делаем вычисления для строк C^T. Я надеялся в конце концов заменить большой цикл цикла вызовом multiprocessing, чтобы распараллелить цикл for, но это будущее. Я также рассмотрю использование Cython для ускорения вычислений.
Вам не нужно mv принимать два аргумента, вместо этого вы можете использовать замыкание (или даже 'functools.partial'): make_mv = lambda c: lambda v: sum (f (c [i] - c [j]) ... '(но вы напишете его с помощью' def', чтобы он был читабельным) и, вероятно, функция должна быть в Cython, или она может быть слишком медленной для такой большой матрицы. – jorgeca
Спасибо! я могу определить mv2 = lambda v: mv (v, C) в моем скрипте, а затем передать mv2 (v) в gmres. Я посмотрю на Cython, так как мой метод довольно медленный для с большими матрицами, с которыми я хочу работать. – user35959
Будьте осторожны, потому что это не закрытие и укусит вас в какой-то момент. Ваша функция 'mv2' будет использовать то, что' C' находится в области видимости, а не та, которая была в области видимости, когда вы создали Если вы хотите, чтобы 'mv2' ссылался на конкретный' C', вам абсолютно необходим один из двух методов выше, и th en 'mv = make_mv (C)'. – jorgeca