Проблема ортогональности в полиномах Лежандра scipy

Недавно я наткнулся на любопытную проблему, касающуюся scipy.special.legendre() (scipy документацию). Полиномы Лежандра должны быть попарно ортогональны. Однако, когда я вычисляю их в диапазоне x=[-1,1] и строю скалярное произведение двух полиномов разной степени, я не всегда получаю нуль или значения, близкие к нулю. Я неправильно истолковываю поведение функций? Далее я написал короткий пример, который производит скалярное произведение некоторых пар полиномов Лежандра:

from __future__ import print_function, division
import numpy as np 
from scipy import special
import matplotlib.pyplot as plt

# create range for evaluation
x = np.linspace(-1,1, 500)

degrees = 6
lp_array = np.empty((degrees, len(x)))

for n in np.arange(degrees):
    LP = special.legendre(n)(x)
    # alternatively:
    # LP = special.eval_legendre(n, x)
    lp_array[n, ] = LP
    plt.plot(x, LP, label=r"$P_{}(x)$".format(n))

plt.grid()
plt.gca().set_ylim([-1.1, 1.1])
plt.legend(fontsize=9, loc="lower right")
plt.show()

График отдельных полиномов на самом деле выглядит нормально: Полиномы Лежандра от 0 до 5 степени

Но если я вычисляю скалярные произведения вручную — умножаю два полинома Лежандра разной степени поэлементно и суммирую их (500 для нормализации)...

for i in range(degrees):
    print("0vs{}: {:+.6e}".format(i, sum(lp_array[0]*lp_array[i])/500))

... Я получаю следующие значения в качестве вывода:

0vs0: +1.000000e+00
0vs1: -5.906386e-17
0vs2: +2.004008e-03
0vs3: -9.903189e-17
0vs4: +2.013360e-03
0vs5: -1.367795e-16

Скалярное произведение первого полинома на себя (как и следовало ожидать) равно единице, а половина других результатов почти равна нулю, но есть некоторые значения порядка 10e-3, и я понятия не имею, почему. Я также попробовал функцию scipy.special.eval_legendre(n, x) - тот же результат :-\

Это ошибка в функции scipy.special.legendre()? Или я что-то не так делаю? Жду конструктивных ответов :-)

привет, Маркус


python scipy polynomials orthogonal
person Markus    schedule 16.09.2016    source источник
comment
спасибо за интеграцию изображения :-)   -  person Markus    schedule 16.09.2016
comment
Хмм хорошо. Я не очень понимаю, почему мой подход не работает правильно, но да, определенно имеет смысл интегрировать полиномиальное произведение. Спасибо за этот момент.   -  person Markus    schedule 16.09.2016
comment
Ваш подход к интеграции не точен. 10^{-3} ~ 1/500 - это порядок величины ошибки, которую вы ожидаете с 500 баллами.   -  person pv.    schedule 16.09.2016
comment
Используйте больше баллов или используйте scipy.integrate.quad. Например. quad(lambda x: special.legendre(0)(x)*special.legendre(2)(x), a=-1, b=1)   -  person Warren Weckesser    schedule 16.09.2016
comment
Что касается моей фактической проблемы (в статистике): здесь у меня есть 126-мерная задача, которая заставляет мои многочлены Лежандра иметь длину 126. Таким образом, это явно слишком мало для scipy.special.legendre, чтобы привести к правильным ортогональным многочленам. Однако это приводит к различным значениям дисперсии при вычислении ковариационной матрицы в зависимости от того, сколько полиномов Лежандра я принимаю во внимание. Есть ли способ справиться с этой проблемой?   -  person Markus    schedule 16.09.2016
comment
Почему полигоны Лежандра должны быть сэмплированы только в 126 точках? Какое отношение количество выборок имеет к размерности задачи? (Если бы к проблеме подошли аналитически, не были бы полигоны функциями от [-1, 1] (бесконечное число точек) до вещественных чисел?)   -  person unutbu    schedule 16.09.2016
comment
Мой набор данных имеет форму (100,60), и поэтому эмпирическая ковариационная матрица имеет форму (100,100). Для уменьшения размеров путем проецирования данных на полиномы Лежандра можно показать с помощью регрессии наименьших квадратов, что оценщик ковариационной матрицы равен S' = P.dot(S_emp).dot(P.T) с P = inv(np.dot(G, G.T)).dot(G.T) и с G, содержащим полиномы Лежандра по столбцам. Для согласованного умножения матриц G должно иметь длину 100, поэтому должны быть и полиномы Лежандра.   -  person Markus    schedule 16.09.2016
comment
К сожалению, у меня нет данных с более высоким разрешением, поскольку они представляют собой среднегодовые значения температуры за период времени в 100 лет.   -  person Markus    schedule 16.09.2016

Ответы (1)


arrow_upward
1
arrow_downward

Как уже отмечали другие, вы получите некоторую ошибку, поскольку выполняете точный интеграл.

Но вы можете уменьшить ошибку, выполняя интеграл как можно лучше. В вашем случае вы все равно можете улучшить свои точки выборки, чтобы сделать интеграл более точным. При выборке используйте «середину» интервалов вместо краев:

x = np.linspace(-1, 1, nx, endpoint=False)
x += 1 / nx   # I'm adding half a sampling interval
              # Equivalent to x += (x[1] - x[0]) / 2

Это дает довольно много улучшений! Если я использую старый метод выборки:

nx = 500
x = np.linspace(-1, 1, nx)

degrees = 7
lp_array = np.empty((degrees, len(x)))

for n in np.arange(degrees):
    LP = special.eval_legendre(n, x)
    lp_array[n, :] = LP

np.set_printoptions(linewidth=120, precision=1)
prod = np.dot(lp_array, lp_array.T) / x.size
print(prod)

Это дает:

[[  1.0e+00  -5.7e-17   2.0e-03  -8.5e-17   2.0e-03  -1.5e-16   2.0e-03]
 [ -5.7e-17   3.3e-01  -4.3e-17   2.0e-03  -1.0e-16   2.0e-03  -1.1e-16]
 [  2.0e-03  -4.3e-17   2.0e-01  -1.3e-16   2.0e-03  -1.0e-16   2.0e-03]
 [ -8.5e-17   2.0e-03  -1.3e-16   1.4e-01  -1.2e-16   2.0e-03  -1.0e-16]
 [  2.0e-03  -1.0e-16   2.0e-03  -1.2e-16   1.1e-01  -9.6e-17   2.0e-03]
 [ -1.5e-16   2.0e-03  -1.0e-16   2.0e-03  -9.6e-17   9.3e-02  -1.1e-16]
 [  2.0e-03  -1.1e-16   2.0e-03  -1.0e-16   2.0e-03  -1.1e-16   7.9e-02]]

Погрешность ~10^-3.

Но, используя «схему выборки средней точки», я получаю:

[[  1.0e+00  -2.8e-17  -2.0e-06  -3.6e-18  -6.7e-06  -8.2e-17  -1.4e-05]
 [ -2.8e-17   3.3e-01  -2.8e-17  -4.7e-06  -2.7e-17  -1.1e-05  -4.1e-17]
 [ -2.0e-06  -2.8e-17   2.0e-01  -5.7e-17  -8.7e-06  -2.3e-17  -1.6e-05]
 [ -3.6e-18  -4.7e-06  -5.7e-17   1.4e-01  -2.1e-17  -1.4e-05  -5.3e-18]
 [ -6.7e-06  -2.7e-17  -8.7e-06  -2.1e-17   1.1e-01   1.1e-17  -2.1e-05]
 [ -8.2e-17  -1.1e-05  -2.3e-17  -1.4e-05   1.1e-17   9.1e-02   7.1e-18]
 [ -1.4e-05  -4.1e-17  -1.6e-05  -5.3e-18  -2.1e-05   7.1e-18   7.7e-02]]

Ошибки теперь ~10^-5 или даже 10^-6, что намного лучше!

person Praveen    schedule 17.09.2016
comment
Во-первых, большое спасибо за ваш ответ. Просто для понимания: добавляя 1/nx, вы просто сдвигаете все тики x на половину интервала. После этого вы на самом деле делаете только скалярное произведение всех полиномов Лежандра друг на друга, верно? Почему тогда вы получаете улучшение условий ошибки? Потому что вы не увеличиваете количество точек выборки. Только полиномы Лежандра оцениваются на немного разных позициях. - person Markus; 18.09.2016