альтернативы или ускорения для инверсии матрицы mpmath

Question

Я пишу некоторый код в python, который требует часто инвертирования больших квадратных матриц (100-200 строк/colums).

Я нахожусь в пределе точности машины, поэтому начал использовать mpmath, чтобы сделать инверсию матрицы произвольной точности, но она очень медленная, даже с использованием gmpy.

Инвертирование случайных матриц размера 20, 30, 60 с точностью 30 (десятичная) занимает ~ 0,19, 0,60 и 4,61 секунды, тогда как те же операции в mathematica составляют 0,0084, 0,015 и 0,055 секунды.

Это использование python3 и mpmath 0.17 (не уверен в версии gmpy) на машине с аркой. Я не уверен, почему mpmath настолько медленнее, но есть ли какая-либо библиотека с открытым исходным кодом, которая будет приближаться к скоростям, которым математика управляет для этого (ровно 1/2, как быстро было бы хорошо)?

Мне не нужна произвольная точность - 128 бит, вероятно, будет достаточно хорошим. Я также просто не понимаю, как mpmath может быть намного медленнее. Он должен использовать совершенно другой алгоритм инверсии матрицы. Чтобы быть конкретным, я использую M**-1.

Есть ли способ заставить его использовать более быстрый алгоритм или ускорить его.

Answer 1

Я предполагаю, что двойная точность не является проблемой для точности конечного результата, но для некоторых матриц это вызывает проблему в промежуточных результатах инверсии. В этом случае давайте рассматривать результат нормальной квантовой (двойной точности) обратной как просто хорошее приближение, а затем использовать это как отправную точку для нескольких итераций метода ньютонов для решения обратной задачи.

Пусть матрица мы пытаемся перевернуть, и X наша оценка обратного. Итерации методы Ньютона просто состоит из:

X = X*(2I - AX) 

Для больших матриц усилия, чтобы вычислить несколько итераций выше почти незначителен по сравнению с отысканием инверсии, и это может значительно повысить точность окончательного результат. Попробуй.

BTW, I - это единичная матрица в приведенном выше уравнении.

EDIT, чтобы добавить код для проверки точности плавающих типов.

Этот код используется для проверки точности поплавкового типа.

x = float128('1.0') 
wun = x 
two = wun + wun 
cnt = 1 
while True: 
    x = x/two 
    y = wun + x 
    if y<=wun: break 
    cnt +=1 

print 'The effective number of mantissa bits is', cnt 

Answer 2

Линера-алгебра в mpmath довольно медленная, к сожалению. Есть много библиотек, которые решают эту проблему намного лучше (например, Sage). Тем не менее, как следует из предложения Стюарта, довольно просто сделать достаточно быстрое высокоточное матричное умножение в Python без установки каких-либо библиотек, используя арифметику с фиксированной точкой. Вот версия с использованием матриц mpmath для ввода и вывода:

def fixmul(A, B, prec): 
    m = A.rows; p = B.rows; n = B.cols; 
    A = [[A[i,j].to_fixed(prec) for j in range(p)] for i in range(m)] 
    B = [[B[i,j].to_fixed(prec) for j in range(n)] for i in range(p)] 
    C = [([0] * n) for r in range(m)] 
    for i in range(m): 
     for j in range(n): 
      s = 0 
      for k in range(p): 
       s += A[i][k] * B[k][j] 
      C[i][j] = s 
    return mp.matrix(C) * mpf(2)**(-2*prec) 

В 256-битной точностью, это умножает две матрицы 200х200 в 16 раз быстрее, чем mpmath для меня. Также нетрудно написать так называемую процедуру инверсии матрицы. Конечно, если элементы матрицы очень большие или очень маленькие, вы хотите сначала их масштабировать.Более твердым решением было бы написать собственные функции матрицы, используя типы с плавающей запятой в gmpy, что должно быть существенно быстрым.

Answer 3

Multiprecision toolbox for MATLAB обеспечивает следующие моменты времени с использованием 128-битной точности (Core i7 930):

20х20 - 0,007 сек

30х30 - 0,019 сек

60x60 - 0,117 сек

200x200 - 3.2 сек

Обратите внимание, что эти цифры намного ниже для современных процессоров.