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

Question

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

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

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

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

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

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

Answer 1

К сожалению, линейная алгебра в 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-битной точностью это умножает две матрицы 200x200 в 16 раз быстрее, чем mpmath для меня. Также нетрудно написать процедуру инверсии матрицы непосредственно таким образом. Конечно, если элементы матрицы очень большие или очень маленькие, вы должны сначала изменить их масштаб. Более твердым решением было бы написать свои собственные матричные функции, используя типы с плавающей точкой в gmpy, что должно быть по существу быстрым.

Answer 2

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

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

X = X*(2I - AX)

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

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

РЕДАКТИРОВАТЬ, чтобы добавить код для проверки точности плавающих типов.

Используйте этот код для проверки точности типа с плавающей точкой.

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 3

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

20x20 - 0,007 с

30x30 - 0,019 с

60x60 - 0,117 сек

200x200 - 3,2 с

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