SciPy: поэлементный неотрицательный метод наименьших квадратов с использованием массива векторов b

Мне нужно решить линейную задачу Ax = b, получив x методом наименьших квадратов. Все элементы x должны быть неотрицательными, поэтому я использую scipy.optimize.nnls (документация здесь).

Проблема в том, что мне нужно много раз решать эту задачу с одной матрицей A и множеством векторов b. У меня есть 3d numpy ndarray, в котором векторы вдоль оси 0 являются векторами b, а две другие оси соответствуют точкам в пространстве. Я хочу вывести все x векторов в соответствующий массив, чтобы сохранить пространственную информацию для каждого ответа.

Первый проход по проблеме выглядит так:

A = np.random.rand(5,3)
b_array = B = np.random.rand(5,100,100)
x_array = np.zeros((3,100,100))

for i in range(100):
    for j in range(100):
        x_array[:,i,j] = sp.optimize.nnls(A, b_array[:,i,j])[0]

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

Некоторое время назад я задал этот очень похожий вопрос о поэлементном умножении матриц . Меня познакомили с np.einsum, который оказался чрезвычайно полезным во многих ситуациях. Я надеялся, что будет аналогичная функция для решения методом наименьших квадратов, но ничего не смог найти. Если кто-нибудь знает о функции, которая может работать, или об альтернативном подходе к эффективному / питоническому решению этой проблемы, это было бы очень признательно!


python numpy scipy arrays least-squares
person Joe    schedule 03.11.2014    source источник

Ответы (1)


arrow_upward
7
arrow_downward

NNLS не имеет решения в закрытой форме, и, кроме совместного использования памяти для матрицы проекта, совместное решение проблем не дает никакого алгоритмического ускорения. Хотя снижение многоцелевой возможности до уровня C может привести к некоторому ускорению, похоже, что реализация scipy поддерживает только одну цель за раз, поэтому зацикливание здесь выглядит как единственный вариант. Проблема смущающе параллельна, поэтому вы можете использовать, например. joblib чтобы распараллелить цикл следующим образом

from joblib import Parallel, delayed
from itertools import product
from scipy.optimize import nnls
results = Parallel(n_jobs=10)(delayed(nnls)(A, b_array[:,i,j])[0]
             for i, j in product(range(100), range(100)))
x_array = np.array(results).reshape(100, 100, -1).transpose(2, 0, 1)

Однако, если вы используете, например. гребневая регрессия или OLS (возможно, бесполезная в вашем случае), тогда решение представляет собой закрытую форму, которую можно получить путем матричного умножения, и все можно сделать за одно изменение формы и матричное умножение, подталкивая многоцелевой аспект вашей проблемы к обработке C-уровня .

person eickenberg    schedule 03.11.2014
comment
Да, я не ожидал алгоритмического ускорения, но надеялся, что будет способ сделать работу C-уровня более эффективной — я плохо понимаю этот аспект языка, поэтому я всегда надеюсь, что есть маленькие хитрости для сделать это быстрее. Я ценю предложение о распараллеливании - мне все равно нужно двигаться в этом направлении, поэтому я попробую. - person Joe; 04.11.2014