Когда не нужно делать собственные разложения для повторного линейного преобразования - PullRequest
Новая
       2

Когда не нужно делать собственные разложения для повторного линейного преобразования

2 голосов
/ 22 января 2020

Допустим, у нас есть точка p, например, (1, 2, 3), к которой мы хотим применить линейное преобразование N раз. Если преобразование обозначается матрицей A, то окончательное преобразование будет иметь вид A^N . p. Умножение матриц дорого, я предполагал, что собственное разложение с последующей диагонализацией ускорит весь процесс. Но, к моему удивлению, этот якобы улучшенный метод занимает больше времени. Что мне здесь не хватает? 

import timeit

mysetup = '''
import numpy as np
from numpy import linalg as LA
from numpy.linalg import matrix_power

EXP = 5     # no. of time linear transformation is applied
LT  = 10    # range from which numbers are picked at random for matrices and points.
N   = 100   # dimension of the vector space

A_init = np.random.randint(LT, size=(N, N))
A = (A_init + A_init.T)/2
p = np.random.randint(LT, size=N)

def run_sim_1():
    An = matrix_power(A, EXP)
    return An @ p

def run_sim_2(): 
    λ, V = LA.eig(A)
    Λ = np.diag(λ)
    Λ[np.diag_indices(N)] = λ ** EXP
    An = V @ Λ @ V.T
    return An @ p
'''

# code snippet whose execution time is to be measured 
# naive implementation
mycode_1 = '''run_sim_1()'''

print(timeit.timeit(setup = mysetup, stmt = mycode_1, number = 1000))
# time taken = 0.14894760597962886

# improved code snippet whose execution time is to be measured
# expecting this to take much less time. 
mycode_2 = '''run_sim_2()'''

# timeit statement 
print(timeit.timeit(setup = mysetup, stmt = mycode_2, number = 1000))
# time taken = 8.035318267997354

Ответы [ 2 ]

2 голосов
/ 22 января 2020

И my_code_1, и my_code_2 содержат только один оператор def. Таким образом, ваши вызовы timeit только определяют время, необходимое для определения функций; функции никогда не называются .

Переместите функцию определения в код настройки и замените операторы, которые должны быть синхронизированы, просто вызовом соответствующей функции, например 

mycode_1 = ''' 
run_sim_1()
'''

Тогда вам следует понизить (на много) значение number, которое вы передаете timeit. И тогда вам нужно будет исправить run_sim_2(), чтобы выполнить правильный расчет: 

def run_sim_2(): 
    λ, V = LA.eig(A)
    Λ = np.diag(λ)
    Λ[np.diag_indices(N)] = λ ** 20
    An = V @ Λ @ V.T
    return An @ p

После внесения этих изменений вы все равно обнаружите, что run_sim_1() быстрее. См. Ответ @ senderle по вероятной причине.

2 голосов
/ 22 января 2020

На этот вопрос сложно ответить авторитетно. Стандартными реализациями умножения матриц и собственного разложения являются O (n ^ 3), поэтому нет априорной причины ожидать, что одно будет быстрее другого. И, как ни странно, мой опыт показывает, что собственное разложение обычно намного медленнее, чем умножение одной матрицы, поэтому этот результат меня не совсем удивляет.

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

# Use binary decomposition to reduce the number of matrix multiplications.
# Here, we iterate over the bits of n, from LSB to MSB, raise `a` to
# increasing powers of 2, and multiply into the result as needed.
z = result = None
while n > 0:
    z = a if z is None else fmatmul(z, z)
    n, bit = divmod(n, 2)
    if bit:
        result = z if result is None else fmatmul(result, z)

Так что на самом деле он не выполняет 20 умножений! Он использует подход «разделяй и властвуй», который уменьшает это число. После продумывания алгоритма, который действительно довольно элегантен, я считаю, что он никогда не будет делать больше, чем 2*log(p) умножений для данной степени p. Этот максимум достигается, когда все биты p равны единице, т. Е. Когда p на единицу меньше степени двойки.

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

Я должен добавить это: не будет ли умножение вектора напрямую быстрее, чем повышение матрицы до сила? Двадцать векторных умножений все равно будут O(n^2), нет? Но, возможно, что вы действительно хотите сделать, это выполнить эту операцию для векторов по 10 Кб, и в этом случае матричный подход к мощности явно превосходит.

...