Матрица Вандермонда

[wataru]

Вас не смущает график числа обусловленности? Странное поведение: линейный рост, а потом какие-то колебания. Выглядят, как ошибка в вычисленях. Возможно, вызванная переполнением переменных.

Заодно, могли бы вы объяснить что вот этот вот код в питоне делает?

n_values[i] = np.array(i)

[kristina_ponomareva]

Здравствуйте!

Нет, не смущает. Колебания вызваны свойствами матрицы Вандермонда, которая часто плохо обусловлена, особенно при больших размерах матрицы. Это приводит к резкому увеличению числа обусловленности и его колебаниям!) ошибок в вычислениях нет.

Что касается вопроса про код, эта строка берет число i, преобразует его в формат массива и добавляется в список n_values. Сейчас, когда вы написали этот комментарий, я поняла, что можно было не делать через array. Но ошибки нет:)

Спасибо за Ваш комментарий!

[wataru]

Колебания вызваны свойствами матрицы Вандермонда

Попробуйте поменять тип данных в матрице на np.float32 (по умолчанию там float64). Посмотрите, что произойдет с вашим графиком. Он станет колебаться раньше.

У вас переполнение. Точности вещественных числел тупо не хватает для вычисления обратной матрицы.

Попробуйте использовать пакет sympy. Или какой-нибудь Pygmp, правда обращение матрицы придется писать руками методом Гаусса, да и вычисление нормы тоже.

Погуглите свойства этой самой матрицы. Уверен, для числа обусловленности есть какая-нибудь аж замкнутая формула.

[Serge3leo]

Попробуйте поменять тип данных в матрице на np.float32 (по умолчанию там float64).

Это да, точность входной матрицы влияет.

Точности вещественных числел тупо не хватает для вычисления обратной матрицы.

А в этом деле, тест с использованием np.float32 практически бесполезен, поскольку NumPy всё равно накапливает скалярное произведение на np.float64.

import numpy as np
import matplotlib.pyplot as plt

def vandermonde_matrix(n, dtype=float, htype=np.longdouble):
    A = np.zeros((n+1,n+1), dtype=dtype)
    for i in range(n+1):
        for j in range(n+1):
            A[i][j] = (htype(i)/htype(n))**htype(j) if n != 0 else 1
    return A

fig, (ax1, ax2) = plt.subplots(nrows=2, sharex=True)

for dt in np.float32, np.float64:
    for ht in np.float32, np.longdouble:
        cond_vandermonde = np.zeros(101, dtype=dt)
        underflow = np.zeros_like(cond_vandermonde)
        rng = range(1, len(cond_vandermonde))
        for i in rng:
            matrix = vandermonde_matrix(i, dtype=dt, htype=ht)
            #cond_vandermonde[i] = np.linalg.cond(matrix)
            U, S, V = np.linalg.svd(matrix)
            cond_vandermonde[i]=np.max(S)/np.min(S)
            #np.testing.assert_approx_equal(np.linalg.cond(U), 1.0)
            #np.testing.assert_approx_equal(np.linalg.cond(V), 1.0)
            #np.testing.assert_approx_equal(np.max(S)/np.min(S),
            #                               np.linalg.cond(np.diag(S)))
            #if cond_vandermonde[i] < 1/np.finfo(dt).eps:
            #    np.testing.assert_approx_equal(np.max(S)/np.min(S),
            #                                   np.linalg.cond(matrix))
            underflow[i] = (i + 1)**2 - i - np.count_nonzero(matrix)
        ax1.semilogy(rng, cond_vandermonde[1:],
                     label=f"m={np.finfo(ht).precision} c={np.finfo(dt).precision}")
        ax2.plot(rng, underflow[1:],
                 label=f"m={np.finfo(ht).precision} c={np.finfo(dt).precision}")

ax1.set_title('Зависимость числа обусловленности от n')
ax1.set_ylabel('Число Обусловленности')
ax1.legend()
ax2.set_ylabel('Число антипереполнений')
ax2.set_xlabel('n')
ax2.legend()
plt.show()

придется писать руками методом Гаусса, да и вычисление нормы тоже.

Я Вас умоляю...

from mpmath import mp
import matplotlib.pyplot as plt

def vandermonde_matrix(n):
    A = mp.zeros(n+1)
    for i in range(n+1):
        for j in range(n+1):
            A[i, j] = (mp.mpf(i)/mp.mpf(n))**j if n else 1
    return A

fig, ax1 = plt.subplots(nrows=1)

for dps in 6, 15, 33:
    mp.dps = dps
    cond_vandermonde = mp.zeros(1, 61)
    rng = range(1, len(cond_vandermonde))
    for i in rng:
        matrix = vandermonde_matrix(i)
        S = mp.svd_r(matrix, compute_uv = False)
        cond_vandermonde[i]=max(S)/min(S)

    ax1.semilogy(rng, cond_vandermonde[1:],
                 label=f"dps={S[0].context.dps}")

ax1.set_title('Зависимость числа обусловленности от n')
ax1.set_ylabel('Число Обусловленности')
ax1.legend()
ax1.set_xlabel('n')
plt.show()

Погуглите свойства этой самой матрицы. Уверен, для числа обусловленности есть какая-нибудь аж замкнутая формула.

Как бы, есть, только она всё равно от сетки зависит. Впрочем, для данной конкретной сетки, на графике и так всё видно.

[Serge3leo]

Колебания вызваны свойствами матрицы Вандермонда

Строго говоря, есть свойства идеальной матрицы Вандермонда и свойства округлённой до заданной точности матрицы Вандермонда, и они отличаются. Отдельный вопрос - сетка, который тоже имеет много гититк, ну пусть останется равномерной:

from mpmath import mp
import matplotlib.pyplot as plt
%config InlineBackend.figure_format='retina'

def vandermonde_matrix(n):
    A = mp.zeros(n+1)
    for i in range(n+1):
        for j in range(n+1):
            A[i, j] = (mp.mpf(i)/mp.mpf(n))**j if n else 1
    return A

fig, ax1 = plt.subplots(nrows=1)

for prec, hiprec in (24, 113), (53, 237), (113, 489):
    mp.prec = prec
    cond_vandermonde = mp.zeros(1, 101)
    cond_vandermonde_hi = mp.zeros(1, 101)
    rng = range(1, len(cond_vandermonde))
    for i in rng:
        mp.prec = prec
        dps = mp.dps
        matrix = vandermonde_matrix(i)
        S = mp.svd_r(matrix, compute_uv = False)
        cond_vandermonde[i]=max(S)/min(S)
        mp.prec = hiprec
        hidps = mp.dps
        S = mp.svd_r(matrix, compute_uv = False)
        cond_vandermonde_hi[i]=max(S)/min(S)

    ax1.semilogy(rng, cond_vandermonde[1:],
                 label=f"vm={dps} cnd={dps}")
    ax1.semilogy(rng, cond_vandermonde_hi[1:],
                 label=f"vm={dps} cnd={hidps}")

ax1.set_title('''Зависимость числа обусловленности от n и от
точности представления матрицы Вандермонда (vm)''')
ax1.set_ylabel('Число обусловленности')
ax1.legend()
ax1.set_xlabel('n')
plt.show()

ошибок в вычислениях нет.

Конечно же, есть. Во-первых, метод вычисления числа обусловленности выбран не самый точный, в NumPy есть быстрая функция np.linalg.cond(matrix), а есть более точное разложение np.linalg.svd(matrix, compute_uv=False).

Легко видеть, что число обусловленности идеальной матрицы Вандермонда неограниченно растёт по одному закону, округленной матрицы Вандермонда тоже растёт неограниченно, но по другому закону точка перегиба: . Ну, а вычисленная оценка числа обусловленности ограничена точностью её вычисления.

[vadimr]

Так ясен пень, что если i делить на 10**17, то результат незначим.

[kristina_ponomareva]

Есть методы интерполяции которые работают всегда, а есть которые нет) про это и рассказываю

[vadimr]

"Всегда" в вещественной арифметике не бывает.

[kristina_ponomareva]

Это да