Estabilidad numérica de polinomios de Zernike de orden superior

Estoy tratando de calcular los momentos de Zernike de orden superior (p. Ej. m=0, n=46) Para alguna imagen. Sin embargo, me encuentro con un problema con respecto al polinomio radial (ver wikipedia ). Este es un polinomio definido en el intervalo [0 1]. Ver el código MATLAB a continuación

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

Sin embargo, esto obviamente tiene problemas numéricos cercanos RHO > 0.9.

Traté de refactorizarlo para polyvalpensar que podría tener algunos mejores algoritmos detrás de escena, pero eso no resolvió nada. La conversión a un cálculo simbólico creó el gráfico deseado, pero fue increíblemente lento incluso para un gráfico simple como se muestra.

¿Existe una forma numéricamente estable de evaluar tales polinomios de alto orden?

Esto se implementa en el siguiente script de 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)

$m = 22$$n = 112$$\rho = 0.7$

$m = 0$$n = 46$1.4e-10

Una posible solución (sugerida por @gammatester) es usar polinomios de Jacobi. Esto evita el problema de la cancelación catastrófica al agregar los grandes coeficientes polinómicos mediante una evaluación polinómica "ingenua".

El polinomio radial de Zernike se puede expresar mediante los polinomios de Jacobi de la siguiente manera (ver ecuación (6) )

Sin embargo, en MATLAB, el uso de jacobiP(n,a,b,x)es inaceptablemente lento para grandes vectores / matrices de x=rho. La jacobiPfunción es en realidad parte de la Caja de herramientas simbólica, y la evaluación del polinomio se difiere al motor simbólico, que intercambia velocidad por precisión arbitraria. Por lo tanto, es necesaria una implementación manual de los polinomios de Jacobi.

$\alpha=m$$\beta=0$$n^*=(n-m/2)$$s$

En MATLAB, esto se traduce en (Jacobi p olice d epartment P olynomial, implementación ' D ' doble ')

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

El polinomio radial real de Zernike es así (para m=abs(m))

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

Nota: esta auto respuesta es solo una solución práctica; siéntase libre de etiquetar otra respuesta que explique por qué esto funciona.