Несколько наборов коэффициентов для наименьших квадратов, соответствующих numpy / scipy

Есть ли способ выполнить несколько одновременных (но не связанных) подгонок методом наименьших квадратов с разными матрицами коэффициентов в numpy.linalg.lstsq или scipy.linalg.lstsq? Например, вот тривиальная линейная аппроксимация, которую я хотел бы сделать с разными значениями x, но с одинаковыми значениями y. В настоящее время мне нужно написать цикл:

x = np.arange(12.0).reshape(4, 3)
y = np.arange(12.0, step=3.0)
m = np.stack((x, np.broadcast_to(1, x.shape)), axis=0)

fit = np.stack(tuple(np.linalg.lstsq(w, y, rcond=-1)[0] for w in m), axis=-1)

Это приводит к набору подгонок с одинаковым наклоном и разными пересечениями, так что fit[n] соответствует коэффициентам m[n].

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

В идеале я хотел бы распространить это на любую «транслируемую» комбинацию a и b, где a.shape[-2] == b.shape[0] точно, а последние измерения должны либо совпадать, либо быть единым (или отсутствовать). Я не особо зацикливаюсь на том, какой размер a представляет различные матрицы: просто было удобно сделать его первым, сократившим цикл.

Есть ли встроенный метод в numpy или scipy, чтобы избежать цикла Python? Меня очень интересует использование lstsq, а не ручное транспонирование, умножение и инвертирование матриц.

Если вы просто хотите избавиться от цикла for, вы можете использовать pinv, а затем матричное умножение: fit=np.linalg.pinv(m)@y

Brenlla 24.10.2018 21:34
Почему в Python есть оператор "pass"?
Почему в Python есть оператор "pass"?
Оператор pass в Python - это простая концепция, которую могут быстро освоить даже новички без опыта программирования.
Некоторые методы, о которых вы не знали, что они существуют в Python
Некоторые методы, о которых вы не знали, что они существуют в Python
Python - самый известный и самый простой в изучении язык в наши дни. Имея широкий спектр применения в области машинного обучения, Data Science,...
Основы Python Часть I
Основы Python Часть I
Вы когда-нибудь задумывались, почему в программах на Python вы видите приведенный ниже код?
LeetCode - 1579. Удаление максимального числа ребер для сохранения полной проходимости графа
LeetCode - 1579. Удаление максимального числа ребер для сохранения полной проходимости графа
Алиса и Боб имеют неориентированный граф из n узлов и трех типов ребер:
Оптимизация кода с помощью тернарного оператора Python
Оптимизация кода с помощью тернарного оператора Python
И последнее, что мы хотели бы показать вам, прежде чем двигаться дальше, это
Советы по эффективной веб-разработке с помощью Python
Советы по эффективной веб-разработке с помощью Python
Как веб-разработчик, Python может стать мощным инструментом для создания эффективных и масштабируемых веб-приложений.
0
1
482
1
Перейти к ответу Данный вопрос помечен как решенный

Ответы 1

Ответ принят как подходящий

Вы можете использовать scipy.sparse.linalg.lsqr вместе с scipy.sparse.block_diag. Я просто не уверен, что это будет быстрее.

Пример:

>>> import numpy as np
>>> from scipy.sparse import block_diag
>>> from scipy.sparse import linalg as sprsla
>>> 
>>> x = np.random.random((3,5,4))
>>> y = np.random.random((3,5))
>>> 
>>> for A, b in zip(x, y):
...     print(np.linalg.lstsq(A, b))
... 
(array([-0.11536962,  0.22575441,  0.03597646,  0.52014899]), array([0.22232195]), 4, array([2.27188101, 0.69355384, 0.63567141, 0.21700743]))
(array([-2.36307163,  2.27693405, -1.85653264,  3.63307554]), array([0.04810252]), 4, array([2.61853881, 0.74251282, 0.38701194, 0.06751288]))
(array([-0.6817038 , -0.02537582,  0.75882223,  0.03190649]), array([0.09892803]), 4, array([2.5094637 , 0.55673403, 0.39252624, 0.18598489]))
>>> 
>>> sprsla.lsqr(block_diag(x), y.ravel())
(array([-0.11536962,  0.22575441,  0.03597646,  0.52014899, -2.36307163,
        2.27693405, -1.85653264,  3.63307554, -0.6817038 , -0.02537582,
        0.75882223,  0.03190649]), 2, 15, 0.6077437777160813, 0.6077437777160813, 6.226368324510392, 106.63227777368986, 1.3277892240815807e-14, 5.36589277249043, array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]))

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

Mad Physicist 24.10.2018 21:53

@MadPhysicist Ага, я просто не знаю, сколько времени теряет алгоритм, потому что ему никто не сказал, что блоки разделены, фактически, что блоки вообще есть. Что ж, я думаю, ты увидишь в своем тайминге.

Paul Panzer 24.10.2018 22:23

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

Mad Physicist 24.10.2018 23:04

Другие вопросы по теме