Каков наиболее точный метод в python для вычисления решения минимальной нормы или решения, полученного из псевдообратного?
Моя цель - решить:
Kc=y
с псевдообратным (то есть <сильным > решением минимальной нормы):
c=K^{+}y
так что модель (надеюсь) - полиномиальная модель высокого уровня f(x) = sum_i c_i x^i. Меня особенно интересует недоопределенный случай, когда мы имеем больше полиномиальных функций, чем данных (несколько уравнений слишком много переменных/неизвестных) columns = deg+1 > N = rows. Примечание K является матрицей вандермодов полиномиальных признаков.
Я изначально использовал функцию python np.linalg.pinv, но потом я заметил, что что-то фанки происходит, как я отметил здесь: Почему разные методы решения Xc = y в python дают разные решения, если они не должны?. В этом вопросе я использую квадратную матрицу для изучения функции на интервале [-1.+1] с полиномом высокой степени. Ответ там предложил мне снизить степень полинома и/или увеличить размер интервала. Основная проблема заключается в том, что неясно, как выбрать интервал или максимальную степень, прежде чем вещь станет ненадежной. Я думаю, что моя основная проблема заключается в том, что выбор такого численно устойчивого диапазона зависит от метода, который я могу использовать. В конце концов, я действительно забочусь о том, что
- метод, который я использую, точно (или действительно близок) к псевдо-обратному для этой задачи полиномиального подбора
- что его численно стабильный
В идеале я хочу попробовать многочлен большой степени, но это может быть ограничено моей машинной точностью. Можно ли увеличить числовую точность машины, используя что-то более точное, чем поплавки?
Кроме того, мне очень важно, что любая функция из используемого python, она обеспечивает самый близкий ответ на псевдо-обратный (и, надеюсь, его численно стабильный, поэтому я действительно могу его использовать). Чтобы проверить, какой ответ для псевдоинверсивного, я написал следующий script:
import numpy as np
from sklearn.preprocessing import PolynomialFeatures
def l2_loss(y,y_):
N = y.shape[0]
return (1/N)*np.linalg.norm(y-y_)
## some parameters
lb,ub = -200,200
N=100
D0=1
degree_mdl = 120
## target function
freq_cos = 2
f_target = lambda x: np.cos(freq_cos*2*np.pi*x)
## evaluate target_f on x_points
X = np.linspace(lb,ub,N) # [N,]
Y = f_target(X) # [N,]
# get pinv solution
poly_feat = PolynomialFeatures(degree=degree_mdl)
Kern = poly_feat.fit_transform( X.reshape(N,D0) ) # low degrees first [1,x,x**2,...]
c_pinv = np.dot(np.linalg.pinv( Kern ), Y)
## get polyfit solution
c_polyfit = np.polyfit(X,Y,degree_mdl)[::-1] # need to reverse to get low degrees first [1,x,x**2,...]
##
c_lstsq,_,_,_ = np.linalg.lstsq(Kern,Y.reshape(N,1))
##
print('lb,ub = {} '.format((lb,ub)))
print('differences with c_pinv')
print( '||c_pinv-c_pinv||^2 = {}'.format( np.linalg.norm(c_pinv-c_pinv) ))
print( '||c_pinv-c_polyfit||^2 = {}'.format( np.linalg.norm(c_pinv-c_polyfit) ))
print( '||c_pinv-c_lstsq||^2 = {}'.format( np.linalg.norm(c_pinv-c_lstsq) ))
##
print('differences with c_polyfit')
print( '||c_polyfit-c_pinv||^2 = {}'.format( np.linalg.norm(c_polyfit-c_pinv) ))
print( '||c_polyfit-c_polyfit||^2 = {}'.format( np.linalg.norm(c_polyfit-c_polyfit) ))
print( '||c_polyfit-c_lstsq||^2 = {}'.format( np.linalg.norm(c_polyfit-c_lstsq) ))
##
print('differences with c_lstsq')
print( '||c_lstsq-c_pinv||^2 = {}'.format( np.linalg.norm(c_lstsq-c_pinv) ))
print( '||c_lstsq-c_polyfit||^2 = {}'.format( np.linalg.norm(c_lstsq-c_polyfit) ))
print( '||c_lstsq-c_lstsq||^2 = {}'.format( np.linalg.norm(c_lstsq-c_lstsq) ))
##
print('Data set errors')
y_polyfit = np.dot(Kern,c_polyfit)
print( 'J_data(c_polyfit) = {}'.format( l2_loss(y_polyfit,Y) ) )
y_pinv = np.dot(Kern,c_pinv)
print( 'J_data(c_pinv) = {}'.format( l2_loss(y_pinv,Y) ) )
y_lstsq = np.dot(Kern,c_lstsq)
print( 'J_data(c_lstsq) = {}'.format( l2_loss(y_lstsq,Y) ) )
используя это, мне удалось заметить, что редко polyfit когда-либо соответствует параметрам, которые использует pinv. Я знаю, что pinv окончательно возвращает псевдоинверсию, поэтому я предполагаю, что если моя главная цель - "убедиться, что я использую псевдоинверсию", то это хорошая идея использовать np.pinv. Тем не менее, я также математически знаю, что псевдоинверс всегда минимизирует ошибку наименьших квадратов J(c) = || Kc - y ||^2 независимо от того, что (доказательство здесь, теорема 11.1.2 446). Таким образом, может быть, моя цель должна состоять в том, чтобы просто использовать функцию python, которая возвращает наименьшую квадратичную ошибку J. Таким образом, я побежал (в недоопределенном случае) сравнение трех методов
- Полигит,
np.polyfit
- псевдо-обратный,
np.linalg.pinv
- наименьшие квадраты,
np.linalg.lstsq
и сравнили, какую ошибку наименьших квадратов они дали мне по данным:
![введите описание изображения здесь]()
Затем я просмотрел странные провалы, которые, как кажется, испытывает функция (что, кстати, похоже на полную тайну, почему есть провалы, если алгоритмы не являются стохастическими), а числа обычно меньше для полифита, например:
lb,ub = (-100, 100)
Data set errors
J_data(c_polyfit) = 5.329753025633029e-12
J_data(c_pinv) = 0.06670557822873546
J_data(c_lstsq) = 0.7479733306782645
учитывая эти результаты, и что псевдоинверс является минимизатором наименьших квадратов, кажется, что лучше всего игнорировать np.pinv. Это лучшее, что нужно сделать? Или я пропущу что-то очевидное?
В качестве дополнительной заметки я вошел в код полифита, чтобы узнать, что именно делает его, имеют более мелкие квадратные ошибки (которые прямо сейчас я использую как способ сказать его наилучшее приближение для псевдо-обратного), и, похоже, он имеет некоторый странный код состояния/числовой стабильности:
# scale lhs to improve condition number and solve
scale = NX.sqrt((lhs*lhs).sum(axis=0))
lhs /= scale
c, resids, rank, s = lstsq(lhs, rhs, rcond)
c = (c.T/scale).T # broadcast scale coefficients
который, как я полагаю, является причиной дополнительной устойчивости к полифиту, который pinv не имеет?
Правильно ли это решение использовать polyfit для моей задачи о приближении линейной системы полиномов высокой степени?
также в этот момент я готов использовать другое программное обеспечение, такое как matlab, если оно дает мне правильную псевдоинверсивную и более численную стабильность (для большинства градусов и любых границ).
Другая случайная идея, которую я только что имел, заключалась в том, что, может быть, есть хороший способ пробовать функцию, так что стабильность псевдо-обратного является хорошей. Я предполагаю, что для аппроксимации косинуса с полиномом требуется некоторый тип числа выборок или расстояние между ними (например, теорема выборки шанквона шеннонов гласит, что базовые функции являются синусоидальными...)
Следует отметить, что, вероятно, инвертирование (или псевдоисследование), а затем умножение - плохая идея. См:
https://www.johndcook.com/blog/2010/01/19/dont-invert-that-matrix/
что речь идет только об обратном, но я предполагаю, что он распространяется и на псевдоинверсы.
прямо сейчас моя путаница заключается в том, что обычно мы не хотим явно вычислять псевдоинверсию явно и делать A^+y=x_min_norm для получения решения минимальной нормы. Однако я бы подумал, что np.lstsq даст ответ, который я хотел, но его ошибка сильно отличается от другой. Я нахожу это чрезвычайно запутанным... заставляю меня думать, что я использую неправильный способ получить решение минимальной нормы в python.
Я не пытаюсь получить регуляризованное решение. Я пытаюсь получить решение минимальной нормы и ничего больше, как можно точнее.
Ответы
Ответ 1
Моя область исследований включает алгоритм сжатия, который по существу называется расширением Фурье. Что является наиболее точным? Он сильно зависит от вектора, который я считаю из-за свойств гладкости. Летом я использовал что-то, называемое Savitsky Golay. Существуют довольно численные и точные способы фильтрации этого. Тем не менее, мой советник имеет метод, который является относительно быстрым и численно стабильным. Область называется продолжением или продолжением Фурье. Как? Я не знаю, разрешено ли мне опубликовать его, вот статья . Если я полагаю, что я, возможно, уже опубликовал летом здесь в python частично.
Это не имеет ничего общего с python, потому что python использует те же базовые библиотеки, что и большинство схем численного кодирования, которые являются BLAS и LAPACK. Netlib находится в режиме онлайн.
Есть несколько других подобных быстрых и численно устойчивых идей, которые могут быть подходящими, я бы рекомендовал. Существует целая книга, посвященная этому b y Boyd. Главы 6 и 7 посвящены этому. Речь идет об общей вариации с регуляризацией из-за основного шума, который у вас может быть в сигнале, который я себе представляю.
Другие аспекты. Возможно, вам придется регулировать СВД из-за плохого кондиционирования. Есть книги, посвященные этому обычно. Просто ответьте на вопрос, какой алгоритм лучше всего. Для алгоритмов существует несколько измерений, и вы не указали на свойства проблемы. Если вы не знали о явление Рунге. Это использование высокопоставленных многочленов является неблагоприятным.
Существует целый класс полиномов Эрмита для решения явлений Гиббса и других методов фильтрации, но это не так хорошо. Вы используете общие функции. Я бы рекомендовал получить Трефтен и Бау. Иногда они делают Революцию Чебычева.
Каково условие числа K. Кроме того, есть что-то, что происходит при подгонке полиномов, называемых явлениями Runges. Вы должны это учитывать. Используйте общие функции, необходимые для приближения низкого ранга для регуляризации, если число условий слишком велико. Я просто прочитал его. Вы используете матрицу Вандермонда. Я собираюсь продемонстрировать это довольно легко. Матрицы Вандермонда плохи. Не используйте их. У них есть узлы.
v = (1:.5:6);
V = vander(v);
c1 = cond(V)
v2 = (1:.5:12);
c2 = cond(vander(v2));
display(c1)
display(c2)
c1 =
6.0469e + 12
c2 =
9.3987e + 32
Я попытался приблизить низкое ранжирование, но матрицы вандермонда не очень приятны. Видеть.
function B = lowapprox(A)
% Takes a matrix A
% Returns a low rank approx of it
% utilizing the SVD
chi = 1e-10;
[U,S,V] = svd(A,0);
DS = diag(S);
aa = find(DS > chi);
s= S(aa,aa);
k = length(aa);
Un = U(:,1:k);
Vn = V(:,1:k)';
B = Un*s*Vn;
end
V2 = vander(v2);
r2 = rank(V2);
c2=cond(V2);
B = lowapprox(V2);
c3 = cond(B);
display(c3)
c2 =
9.3987e+32
c3 =
3.7837e+32
ничего не делает... Если вы не знаете, что происходит, когда вы получаете это обратное, число условий равно максимальному сингулярному значению над минимумом, так что у вас есть очень маленькие особые значения при точности машины.
Кроме того, я думаю, у вас есть некоторая путаница в отношении минимальной нормы и регуляризации. Вы сказали, что хотите минимальную норму в смысле наименьших квадратов. SVD дает наименьшие квадраты.
Свойство 9, A построено из базы SVD. Это покрыто трефетеном, но матрица вандермонда плохо обусловлена.
даже небольшие плохо построенные матрицы вандермонда будут терять его. Теперь о вашем решении. Не используйте матрицы вандермонда. Построим полином иначе. Лучшая идея - барицентрическая интерполяция лагранжа. Библиотека здесь
Вот пример в Matlab.
t= (0:.01:pi)';
f = cos(t);
data = [t,f];
f1 = barylag(data,t)
display(err =norm(f1-f1))
err =
0
barylag берется с сайта Matlab. Поскольку я не могу действительно прокомментировать ваши несоответствия, вы должны понять, как именно выполняется lsqr. Алгоритмы Lsqr - это методы Крылова. Это описано в Трефетене. SVD также У меня есть пример на моей странице quora о числовой стабильности с QR, который является тем, как вы на самом деле строите эти алгоритмы
