Stabilność numeryczna wielomianów zerowych wyższego rzędu

Próbuję obliczyć wyższego rzędu (np m=0, n=46) chwile Zernike jakiegoś obrazu. Mam jednak problem dotyczący wielomianu promieniowego (patrz wikipedia ). Jest to wielomian zdefiniowany w przedziale [0 1]. Zobacz kod MATLAB poniżej

function R = radial_polynomial(m,n,RHO)
    R = 0;
    for k = 0:((n-m)/2)        
        R = R + (-1).^k.*factorial(n-k) ...
            ./ ( factorial(k).*factorial((n+m)./2-k) .* factorial((n-m)./2-k) ) ...
            .*RHO.^(n-2.*k);
    end
end

To oczywiście napotyka na problemy numeryczne w pobliżu RHO > 0.9.

Próbowałem przeredagować go, tak aby polyvalmyślał, że może mieć lepsze algorytmy za kulisami, ale to niczego nie rozwiązało. Przekształcenie go w obliczenie symboliczne stworzyło pożądany wykres, ale było zadziwiająco powolne, nawet w przypadku prostego wykresu, takiego jak pokazano.

Czy istnieje stabilny numerycznie sposób oceny wielomianów wysokiego rzędu?

matlab polynomials numerical-limitations

— Sanchises
źródło

Często lepiej jest stosować wielomiany ortogonalne, tutaj wielomiany Jacobiego . Próbowałeś już mathworks.com/help/symbolic/jacobip.html i relacji

R_{n}^{m} (r) = (- 1)^{(n - m) / 2} r^{m} P_{(n - m) / 2}^{(m, 0)} (1 - 2 r^{2}) ?

$R_n^{\,m}(r) = (-1)^{(n-m)/2}\,r^m\,P_{(n-m)/2}^{\,(m,0)}(1-2r^2)?$

— gammatester

@gammatester To działa! Czy mógłby Pan wyjaśnić, dlaczego tak się dzieje?

— Sanchises

Miło słyszeć, że to działa. Niestety nie mogę udzielić zdecydowanej odpowiedzi z dwóch powodów. Po pierwsze: chociaż powszechnie wiadomo, że wielomiany ortogonalne mają lepsze właściwości stabilności niż forma standardowa, nie znam formalnego dowodu (szczególnie w tym przypadku). Po drugie, nie używam Matlaba i nie mogę podać danych dla zaimplementowanych wielomianów Jacobiego.

— gammatester

@ Sanchises Nie ma tu darmowego lunchu: tylko dlatego, że coś jest wielomianem, nie oznacza, że bezpośrednia formuła pod względem mocy jest właściwym sposobem na jego obliczenie, a dokładne obliczenie wielomianów Jacobiego nie jest samo w sobie trywialne - nie robisz tego przez współczynniki, więc nie jest tak tanie.

— Kirill,

Powodem użycia wielomianów Jacobi jest to, że pozbyłeś się katastrofalnego anulowania w swojej formule (spójrz na wszystkie te czynniki oscylacyjne o bardzo dużych współczynnikach!), A domyślna procedura oceny wielomianu Jacobi jest implementowana ostrożnie w bibliotece, więc jest gwarantowana być punktualnym. Większość pracy tutaj polega na upewnieniu się, że wielomiany Jacobi są dokładnie ocenione.

— Kirill

Odpowiedzi:

W tym artykule Honarvar i Paramesran opracowali interesującą metodę obliczenia wielomianów promieniowych Zernike w bardzo przyjemny sposób rekurencyjny. Formuła rekurencyjna jest zaskakująco prosta, bez dzielenia lub mnożenia przez duże liczby całkowite:

R_{n}^{m} (ρ) = ρ (R_{n - 1}^{| m - 1 |} (ρ) + R_{n - 1}^{m + 1} (ρ)) - R_{n - 2}^{m} (ρ)

$R^m_n(\rho) = \rho \left(R^{|m-1|}_{n-1}(\rho)+R^{m+1}_{n-1}(\rho)\right) - R^{m}_{n-2}(\rho)$ Poleciłbym rzucić okiem na rysunek 1 w pracy Honarvara i Paramesrana, który wyraźnie ilustruje zależności między różnymi wielomianami Zernike.

Jest to zaimplementowane w następującym skrypcie Octave:

clear                                     % Tested with Octave instead of Matlab
N = 120;
n_r = 1000;
R = cell(N+1,N+1);
rho = [0:n_r]/n_r;
rho_x_2 = 2*[0:n_r]/n_r;

R{0+1,0+1} = ones(1,n_r+1);               % R^0_0  Unfortunately zero based cell indexing is not possible
R{1+1,1+1} = R{0+1,0+1}.*rho;             % R^1_1  ==>  R{...+1,...+1} etc.
for n = 2:N,
    if bitget(n,1) == 0,                  % n is even
        R{0+1,n+1} = -R{0+1,n-2+1}+rho_x_2.*R{1+1,n-1+1};                % R^0_n
        m_lo = 2;
        m_hi = n-2;
    else
        m_lo = 1;
        m_hi = n-1;
    end
    for m = m_lo:2:m_hi,
        R{m+1,n+1} = rho.*(R{m-1+1,n-1+1}+R{m+1+1,n-1+1})-R{m+1,n-2+1};  % R^m_n
    end
    R{n+1,n+1} = rho.*R{n-1+1,n-1+1};                                    % R^n_n
end;


Z = @(m,n,rho) (-1)^((n-m)/2) * rho.^m .* jacobiPD((n-m)/2,m,0,1-2*rho.^2);
m = 22;
n = 112;
figure
plot(rho,Z(m,n,rho))
hold on
plot(rho,R{m+1,n+1},'r');
xlabel("rho")
ylabel("R^{22}_{112}(rho)")
legend("via Jacobi","recursive");
%print -djpg plt.jpg

m = 0;
n = 46;
max_diff_m_0_n_46 = norm(Z(m,n,rho)-R{m+1,n+1},inf)

Na przykład liczba wygenerowana przez ten kod pokazuje, że za pomocą $m = 22$ , i $n = 112$ , katastrofalne anulowanie następuje w pobliżu $\rho = 0.7$ , jeśli wielomiany promieniowe Zernike są obliczane za pomocą wielomianów Jacobiego. Dlatego należy również martwić się o dokładność wielomianów Zernike niższego stopnia.

Metoda rekurencyjna wydaje się znacznie bardziej odpowiednia do stabilnego obliczania wielomianów zerowych wyższego rzędu. Niemniej jednak dla $m = 0$ i $n = 46$ , maksymalna różnica między metodą Jacobi a metodą rekurencyjną wynosi (tylko?) 1.4e-10, co może być wystarczająco dokładne dla Twojej aplikacji.

— wim
źródło

Twoja fabuła wygląda jak błąd w Matlabie jacobiPD, a nie jak jakakolwiek katastroficzna anulacja.

— Kirill,

@Kiril: Użyłem Sanchises JacobiPDz jego odpowiedzi . Działa to dobrze w przypadku wielomianów niskiego rzędu. Na przykład za pomocą

n = 30

$n=30$ , arbitralne

m

$m$ i arbitralne

ρ

$\rho$ , różnica między dwiema metodami jest mniejsza niż 6.9e-13. Chociaż poszczególne warunki w podsumowaniu JacobiPDsą małe, mogą się zwiększyć po pomnożeniu przez factorial(n+a) * factorial(n+b). Ponadto mają naprzemienny znak, który jest idealnym przepisem na katastrofalne anulowanie.

— wim,

(ciąg dalszy) Np. z

m = 22

$m=22$ i

n = 112

$n=112$ wyrażenie

1/(factorial(s)*factorial(n+a-s)*factorial(b+s)*factorial(n-s)) *              ((x-1)/2).^(n-s).*((x+1)/2).^s * factorial(n+a) * factorial(n+b)

może stać się tak duże, jak 1.4e18, a suma jest dopiero w -2.1końcu. Możesz nazwać to błędem, ale z nieskończoną precyzją odpowiedź byłaby poprawna. Czy możesz wyjaśnić, co rozumiesz przez „brak ogólnej katastrofalnej rezygnacji”?

— wim,

@ wim Nie zauważyłem, że to nie Matlab. Jeśli czyjaś wielomianowa implementacja Jacobi jest wystarczająco dobra do ich celu, nie ma problemu. Powiedziałem tylko, że to błąd, ponieważ źle zrozumiałem i pomyślałem, że jest to funkcja wbudowana (spodziewam się, że funkcje biblioteki będą bardzo solidne). Przez „ogólny” miałem na myśli, że jeśli nie wiesz, jak ta funkcja jest zaimplementowana, nie możesz nazwać niepoprawnych danych wyjściowych „katastroficznym anulowaniem”, jak ogólny termin dla wszystkich rodzajów błędów, ale to tylko moje nieporozumienie kod działał.

— Kirill,

Żeby było jasne: mój kod nie jest rekurencyjny. Jest to iteracyjna standardowa trzyterminowa relacja powtarzalności (podobna do wielomianów Czebyszewa), która ma być normalnie bardziej stabilna niż np. Postać Hornera dla wielomianów.

— gammatester

Możliwym rozwiązaniem (sugerowanym przez @gammatester) jest użycie wielomianów Jacobi. To omija problem katastroficznego anulowania przy dodawaniu dużych współczynników wielomianowych poprzez „naiwną” ocenę wielomianu.

Promieniowy wielomian Zernike może być wyrażony przez wielomiany Jacobiego w następujący sposób (patrz równanie (6) )

R_{n}^{m} (ρ) = (- 1)^{(n - m) / 2} ρ^{m} \cdot P_{(n - m) / 2}^{(m, 0)} (1 - 2 ρ^{2})

$R^m_n(\rho) = (-1)^{(n-m)/2}\rho^m \cdot P^{(m,0)}_{(n-m)/2} \Big(1-2\rho^2 \Big)$

Jednak w MATLAB użycie jacobiP(n,a,b,x)jest niedopuszczalnie powolne w przypadku dużych wektorów / macierzy x=rho. Ta jacobiPfunkcja jest właściwie częścią Symbolicznego zestawu narzędzi, a ocena wielomianu jest odraczana do silnika symbolicznego, który zamienia prędkość na dowolną precyzję. Konieczna jest zatem ręczna implementacja wielomianów Jacobi.

Ponieważ wszystkie parametry funkcji Jacobi są nieujemne ( $\alpha=m$ , $\beta=0$ , $n^*=(n-m/2)$ ), możemy użyć następującego wyrażenia (patrz Wikipedia , zauważ, że wpisałem wartości dla $s$ )

\begin{matrix} P_{n}^{(α, β)} (ρ) = (n + α)! (n + β)! \cdot \\ \sum_{s = 0}^{n} [\frac{1}{s! (n + α - s)! (β + s)! (n - s)!} {(\frac{x - 1}{2})}^{n - s} {(\frac{x + 1}{2})}^{s}] \end{matrix}

$\begin{multline} P_n^{(\alpha,\beta)}(\rho) = (n+\alpha)!(n+\beta)! \quad \cdot\\ \sum_{s={0}}^n \left[\frac{1}{s!(n+\alpha-s)!(\beta+s)!(n-s)!} \left(\frac{x-1}{2}\right)^{n-s}\left(\frac{x+1}{2}\right)^{s} \right]\end{multline}$

MATLAB przekłada się (Jacobim ~~P olice d epartment~~ P olynomial ' D ouble' Wykonanie)

function P = jacobiPD(n,a,b,x)
    P = 0;
    for  s  0:n
        P = P + ...
            1/(factorial(s)*factorial(n+a-s)*factorial(b+s)*factorial(n-s)) * ...
            ((x-1)/2).^(n-s).*((x+1)/2).^s;
    end
    P = P*factorial(n+a) * factorial(n+b);
end

Rzeczywisty wielomian promieniowy Zernike jest zatem (dla m=abs(m))

Z = @(m,n,rho) (-1)^((n-m)/2) * rho.^m .* jacobiPD((n-m)/2,m,0,1-2*rho.^2);

Uwaga: ta odpowiedź jest tylko praktycznym rozwiązaniem; dodaj inną odpowiedź, która wyjaśnia, dlaczego to działa.

— Sanchises
źródło