c++ - tercer - Adaptar datos a un polinomio de 3er grado

polinomio de grado 3 ejemplo (3)

Actualmente estoy escribiendo un programa en C ++ donde tengo vectores de datos independientes y dependientes que me gustaría adaptar a una función cúbica. Sin embargo, tengo problemas para generar un polinomio que se ajuste a mis datos.

Parte del problema es que no puedo usar varios paquetes numéricos, como GSL (historia larga); es posible que incluso sea excesivo para mi caso. No necesito una solución muy general para el ajuste de mínimos cuadrados. Específicamente, quiero ajustar mis datos a una función cúbica. Tengo acceso a la biblioteca de vectores de Sony, que admite matrices 4x4 y puede calcular sus inversas, entre otras cosas.

Mientras prototipé esto en Scilab, utilicé una función como:

function p = polyfit(x, y, n) m = length(x); aa = zeros(m, n+1) aa(:,1) = ones(m,1) for k = 2:n+1 aa(:,k) = x.^(k-1) end p = aa/y endfunction

Desafortunadamente, esto no se adapta bien a mi entorno actual. El ejemplo anterior necesita soportar una matriz de dimensiones M x N + 1. En mi caso, eso es M x 4, donde M depende de la cantidad de datos de muestra que tengo. También está el problema de la división izquierda. Necesitaría una biblioteca de matriz que admite el inverso de matrices de dimensiones arbitrarias.

¿Hay un algoritmo para mínimos cuadrados donde puedo evitar tener que calcular aa / y, o al menos limitarlo a una matriz de 4x4? Supongo que estoy tratando de reescribir el algoritmo anterior en un caso más simple que funcione para ajustarse a un polinomio cúbico. No estoy buscando una solución de código, pero si alguien puede señalarme en la dirección correcta, lo agradecería.

Nunca debe hacer una inversión de matriz completa por razones de estabilidad. Hacer la descomposición LU y la sustitución hacia adelante y hacia atrás. Las otras soluciones son acertadas de otra manera.

Sí, podemos limitar el problema a la informática con "una matriz 4x4". El ajuste por mínimos cuadrados de un punto cúbico, incluso para M puntos de datos, solo requiere la solución de cuatro ecuaciones lineales en cuatro incógnitas. Suponiendo que todas las coordenadas x son distintas, la matriz de coeficientes es invertible, por lo que, en principio, el sistema puede resolverse invirtiendo la matriz de coeficientes. Suponemos que M es más de 4, como sería típicamente el caso para los ajustes de mínimos cuadrados.

Aquí hay un resumen de Maple, Cómo ajustar un cubículo a los datos , que oculta casi por completo los detalles de lo que se está resolviendo. Los criterios mínimos de primer orden (primeras derivadas con respecto a los coeficientes como parámetros del error de suma de cuadrados) nos dan las cuatro ecuaciones lineales, a menudo denominadas ecuaciones normales .

Puede "ensamblar" estas cuatro ecuaciones en código, luego aplicar su matriz inversa o una estrategia de solución más sofisticada. Obviamente, necesita tener los puntos de datos almacenados de alguna forma. Una posibilidad es dos matrices lineales, una para las coordenadas xy una para las coordenadas y, ambas de longitud M el número de puntos de datos.

NB: Voy a discutir este conjunto de matriz en términos de subíndices de matriz basados en 1. Los coeficientes polinomiales son en realidad una aplicación donde los subíndices de matriz basados en 0 hacen las cosas más simples y limpias, pero reescribirla en C o en cualquier otro idioma que favorezca los subíndices basados en 0 se deja como un ejercicio para el lector.

El sistema lineal de ecuaciones normales se expresa más fácilmente en forma de matriz refiriéndose a una matriz A de Mx4 cuyas entradas son potencias de datos de coordenadas x:

A (i, j) = coordenada x del ith punto de datos elevado a la potencia j-1

Deje que A ''denote la transposición de A, de modo que A''A sea una matriz de 4x4.

Si dejamos d ser una columna de longitud M que contiene las coordenadas y de los puntos de datos (en el orden dado), entonces el sistema de ecuaciones normales es solo esto:

A''A u = A ''d

donde u = [p0, p1, p2, p3] ''es la columna de coeficientes para el polinomio cúbico con ajuste de mínimos cuadrados:

P (x) = p0 + p1 * x + p2 * x ^ 2 + p3 * x ^ 3

Sus objeciones parecen centrarse en una dificultad para almacenar y / o manipular la matriz A de Mx4 o su transposición. Por lo tanto, mi respuesta se centrará en cómo ensamblar la matriz A''A y la columna A''d sin almacenar explícitamente todas las A a la vez. En otras palabras, haremos las multiplicaciones indicadas matriz-matriz y matriz-vector implícitamente para obtener un sistema 4x4 que pueda resolver:

C u = f

Si piensas que la entrada C (i, j) es el producto de la i-ésima fila de A ''con la j-ésima columna de A, más el hecho de que la i-ésima fila de A'' es realmente la transposición de la i-ésima columna de A , debe quedar claro que:

C (i, j) = SUM x ^ (i + j-2) sobre todos los puntos de datos

¡Este es ciertamente un lugar donde la exposición se simplificaría al usar subíndices basados en 0!

Puede tener sentido acumular las entradas para la matriz C, que depende solo del valor de i + j, es decir, una llamada matriz de Hankel , en una matriz lineal de longitud 7 tal que:

W (k) = SUM x ^ k sobre todos los puntos de datos

donde k = 0, .., 6. La matriz 4x4 C tiene una estructura "a rayas" que significa que solo aparecen estos siete valores. Al pasar por la lista de coordenadas x de puntos de datos, puede acumular las contribuciones apropiadas de cada potencia de cada punto de datos en la entrada correspondiente de W.

Se puede usar una estrategia similar para ensamblar la columna f = A ''d, es decir, recorrer los puntos de datos y acumular las siguientes cuatro sumas:

f (k) = SUMA (x ^ k) * y sobre todos los puntos de datos

donde k = 0,1,2,3. [Por supuesto, en las sumas anteriores, los valores x, y son las coordenadas de un punto de datos común.]

Advertencias: Esto satisface el objetivo de trabajar solo con una matriz 4x4. Sin embargo, uno normalmente trata de evitar la formación explícita de la matriz de coeficientes para las ecuaciones normales porque estas matrices son a menudo lo que en el análisis numérico se llama mal condicionado. En particular, los casos en que las coordenadas x están muy cercadas pueden causar dificultades cuando uno trata de resolver el sistema invirtiendo la matriz de coeficientes.

Un enfoque más sofisticado para resolver estas ecuaciones normales sería el método de gradiente conjugado en las ecuaciones normales , que se puede hacer con un código que compute los productos de matriz vectorial A u y A ''v una entrada a la vez (usando lo que dijimos arriba) sobre las entradas de A).

La precisión del método de gradiente conjugado a menudo es satisfactoria debido a su enfoque iterativo natural, esp. cuando uno puede calcular los productos de puntos requeridos con un poco de precisión extra.

Aquí está la página en la que estoy trabajando, aunque esa página en sí misma no aborda su pregunta directamente. El resumen de mi respuesta sería:

Si no puede trabajar con matrices Nx4 directamente, haga esos cálculos matriciales "manualmente" hasta que tenga el problema reducido a algo que tenga solo matrices 4x4 o más pequeñas. En esta respuesta, explicaré cómo hacer los cálculos de matriz específicos que necesita "manualmente".

Supongamos que tiene un montón de puntos de datos (x1,y1)...(xn,yn) y está buscando la ecuación cúbica y = ax^3 + bx^2 + cx + d que mejor se ajuste a esos puntos.

Luego, siguiendo el enlace de arriba, escribirías esta ecuación:

Escribiré A , x y B para esas matrices. Luego, siguiendo mi enlace de arriba, te gustaría multiplicar por la transposición de A , que te dará la matriz de 4x4 A ^T * A que puedes invertir. En ecuaciones, el siguiente es el plan:

A * x = B .................... [con lo que comenzamos]
(A ^T * A) * x = A ^T * B ..... [multiplicar por A ^T ]
x = (A ^T * A) ^-1 * A ^T * B ... [multiplicar por el inverso de A ^T * A]

Dijiste que estás contento con la inversión de matrices 4x4 , así que si podemos codificar una forma de llegar a estas matrices sin utilizar realmente objetos de matriz, deberíamos estar bien.

Entonces, este es un método, aunque podría ser demasiado parecido a crear su propia biblioteca de matriz para su gusto. :)

Escribe una ecuación explícita para cada una de las 16 entradas de la matriz 4x4. La entrada (i,j) th (empiezo con (0,0)) viene dada por x ₁ ⁱ * x ₁ ^j + x ₂ ⁱ * x ₂ ^j + ... + x _N ⁱ * x _N ^j .
Invierta esa matriz 4x4 usando su biblioteca de matriz. Eso es ( A ^T * A ) ^-1 .
Ahora todo lo que necesitamos es A ^T * B , que es una matriz 4x1. La i-ésima entrada está dada por x ₁ ⁱ * y ₁ + x ₂ ⁱ * y ₂ + ... + x _N ⁱ * y _N.
Multiplica nuestra matriz 4x4 creada a mano ( A ^T * A ) ^-1 por nuestra matriz 4x1 creada a mano A ^T * B para obtener la matriz 4x1 de coeficientes de mínimos cuadrados para tu cubo.

¡Buena suerte!