Numpy / Scipy: решение нескольких наименьших квадратов с одинаковой матрицей - PullRequest
Новая
       65

Numpy / Scipy: решение нескольких наименьших квадратов с одинаковой матрицей

1 голос
/ 14 февраля 2020

Я сталкиваюсь с проблемой наименьших квадратов, которую решаю через scipy.linalg.lstsq(M,b), где:

  • M имеет форму (n,n)
  • b имеет форму (n,)

Проблема в том, что я должен решить ее кучу времени для разных b. Как я могу сделать что-то более эффективное? Я предполагаю, что lstsq делает много вещей независимо от значения b.

Идеи?

Ответы [ 3 ]

1 голос
/ 15 февраля 2020

В случае, если ваша линейная система хорошо определена, я сохраню M LU-разложение и буду использовать его для всех b по отдельности или просто сделаю один вызов решения для 2d-массива B, представляющего горизонтально расположенные b, это действительно зависит от вашей проблемы, но это глобально та же идея. Предположим, у вас есть каждый b по одному, затем: 

import numpy as np
from scipy.linalg import lstsq, lu_factor, lu_solve, svd, pinv

# as you didn't specified any practical dimensions
n = 100
# number of b's
nb_b = 10

# generate random n-square matrix M
M = np.random.rand(n**2).reshape(n,n)

# Set of nb_b of right hand side vector b as columns
B = np.random.rand(n*nb_b).reshape(n,nb_b)

# compute pivoted LU decomposition of M
M_LU = lu_factor(M)
# then solve for each b
X_LU = np.asarray([lu_solve(M_LU,B[:,i]) for i in range(nb_b)])

, но если он недооценен или переопределен, вам нужно использовать lstsq , как вы это сделали : 

X_lstsq = np.asarray([lstsq(M,B[:,i])[0] for i in range(nb_b)])

или просто сохранить псевдообратную M_pinv с pinv (построена на lstsq) или pinv2 (построена на SVD): 

# compute the pseudo-inverse of M
M_pinv = pinv(M)
X_pinv = np.asarray([np.dot(M_pinv,B[:,i]) for i in range(nb_b)])

или вы также можете выполнить работу самостоятельно, как, например, в pinv2, просто сохраните SVD в M и решите это вручную: 

# compute svd of M
U,s,Vh = svd(M)

def solve_svd(U,s,Vh,b):
    # U diag(s) Vh x = b <=> diag(s) Vh x = U.T b = c
    c = np.dot(U.T,b)
    # diag(s) Vh x = c <=> Vh x = diag(1/s) c = w (trivial inversion of a diagonal matrix)
    w = np.dot(np.diag(1/s),c)
    # Vh x = w <=> x = Vh.H w (where .H stands for hermitian = conjugate transpose)
    x = np.dot(Vh.conj().T,w)
    return x

X_svd = np.asarray([solve_svd(U,s,Vh,B[:,i]) for i in range(nb_b)])

, которые все дают тот же результат, если проверено с помощью np.allclose (если система не является точно определенной, что приводит к отказу прямого подхода LU). Наконец, с точки зрения производительности: 

%timeit M_LU = lu_factor(M); X_LU = np.asarray([lu_solve(M_LU,B[:,i]) for i in range(nb_b)])
1000 loops, best of 3: 1.01 ms per loop

%timeit X_lstsq = np.asarray([lstsq(M,B[:,i])[0] for i in range(nb_b)])
10 loops, best of 3: 47.8 ms per loop

%timeit M_pinv = pinv(M); X_pinv = np.asarray([np.dot(M_pinv,B[:,i]) for i in range(nb_b)])
100 loops, best of 3: 8.64 ms per loop

%timeit U,s,Vh = svd(M); X_svd = np.asarray([solve_svd(U,s,Vh,B[:,i]) for i in range(nb_b)])
100 loops, best of 3: 5.68 ms per loop

Тем не менее, вы должны проверить их с соответствующими размерами.

Надеюсь, это поможет.

0 голосов
/ 15 февраля 2020

Вы можете разделить M на продукты QR или SVD и найти решение lsq вручную.

0 голосов
/ 14 февраля 2020

Ваш вопрос неясен, но я предполагаю, что вы хотите вычислить уравнение от Mx=b до scipy.linalg.lstsq(M,b) для различных массивов (b0, b1, b2..). Если это так, вы можете просто распараллелить процесс с concurrent.futures.ProcessPoolExecutor. Документация для этого довольно проста и может помочь python запустить несколько решателей scipy одновременно.

Надеюсь, это поможет.

...