Recibido: 13-06-2021 / Revisado: 22-06-202=
1 /
Aceptado: 11-07-2021 / Publicado: 05-08-2021
Usando
la Factorización de Rango Completo de Cholesky en la Solución de los Mínimos
Cuadrados Ponderados
DOI:
Full Rank Cholesky Factorization for the Re-Weighted Least Squares Applications<= o:p>
Zenaida Natividad Castillo Marrero. [1], Ernesto Antonio Ponsot
Balaguer. [2],
Franklin José Camacho. [3] & María Victoria León Sánchez. [4]
Factorization of matrices are required
frequently in sciences research, in order to simplify computations when sol=
ving
linear or non-linear systems of equations. In particular, these systems ari=
se in
data analysis and also after modeling industrial process, playing an import=
ant
role in the final solution. When the matrices involved are full rank and sq=
uared
classical factorizations can be used efficiently. However, when one has to =
deal
with rectangular and/or rank deficient matrices, extensions of these
factorizations should be considered, if it is possible with the same good p=
roperties.
Objective: In this work a full rank Cholesky factorization for rank
deficient matrices is described, along with its use in a particular applica=
tion
arising in Statistics, where usually one has to solve overdetermined system=
s of
equations via least squared method or its extensions. Methodology: We
describe an algorithm to generate this factorization in an iteratively
reweighed least squared method applied to maximum likelihood parameters
estimation. Results: Preliminary experiments are presented over
synthetic data, supporting the use of this decomposition for rectangular or
rank deficient matrices arising in similar applications. Conclusions=
: The
proposal of using the full rank Cholesky factorization in the iteratively
reweighted least squared to estimate maximum likehood parameters is promisi=
ng
and could be used to improve CPU time in the solution of similar problems. =
Resumen
Es frecuente encontrar en la investigación en las ciencias e
ingenierías, que se requiere de factorizaciones de matrices para facilitar o
simplificar los cálculos, ya sea para resolver sistemas de ecuaciones linea=
les,
o no lineales, provenientes de análisis de datos o de modelado o simulación=
de
procesos industriales. Cuando las matrices son cuadradas y de rango completo
existen factorizaciones que permiten hacer este trabajo eficientemente; sin
embargo, cuando las matrices son rectangulares y/o deficientes en rango, las
extensiones de estas factorizaciones se hacen necesarias y deben ser tratad=
as
cuidadosamente para que los resultados sean confiables. Objetivo: En este trabajo se describe la descomposición de rango completo de
Cholesky y su uso en algunas aplicaciones en áreas como la Estadística, don=
de a
menudo se requiere resolver sistemas sobredeterminados por el método de los
mínimos cuadrados o extensiones del mismo. Metodología. Se describen=
los
pasos para generar esta descomposición y su uso en la estimación máximo
verosímil de parámetros por el método de los mínimos cuadrados reponderados=
. Resultados:
Se presenta experimentación preliminar en datos sintéticos, con resultados =
que
sustentan el uso de esta descomposición para matrices rectangulares o defic=
ientes
en rango. =
Conclusión:
La factorización propuesta para resolver el problema de estimación de
parámetros con la función de verosimilitud es satisfactoria y pudiera ser
utilizada para optimizar el tiempo de cómputo de la solución de problemas s=
imilares.
Palabras claves: Factori=
zación
de Matrices, Deficiencia en rango, Factorización de rango completo, Cholesk=
y,
Mínimos Cuadrados Ponderados.
Uno de los cálculos que aparece con mayor frecuencia en la
resolución de problemas de modelado industrial, o en el análisis de datos, =
es
sin duda la solución de sistemas lineales de la forma
Cuando la matriz es simétrica y definida positiva es convenien=
te
usar la descomposición de Cholesky, ya que se disminuye considerablemente el
tiempo de cálculo y se aprovecha la fortaleza de esta factorización. En la
mayoría de los lenguajes de programación que se manejan en la actualidad, e=
xisten
librerías con códigos robustos para estas descomposiciones, los cuales func=
ionan
apropiadamente para aplicaciones generales densas o dispersas, ver Datta (2=
010)
para mayor referencia.
Ahora bien, en ocasiones las matrices tienen una estructura
particular que pudiera explotarse para ganar precisión o rapidez en los
cálculos; en este caso los investigadores prefieren elaborar sus propios
códigos.
Para matrices semidefinidas positivas y para matrices
rectangulares de rango incompleto aún es posible hallar una extensión de la
factorización de Cholesky con las mismas propiedades o el mismo potencial. =
En este trabajo se describen este tipo d=
e factorizaciones
y se presenta su uso en aplicaciones de interés en áreas de proyección en la
actualidad. La experimentación preliminar muestra la efectividad de la
propuesta en un ejemplo de aplicación.
En esta sección y para efectos de comparación se presenta un
resumen de las factorizaciones clásicas de matrices, que han resultado de m=
ucha
utilidad para resolver sistemas lineales de tamańo moderado cuya matriz de
coeficientes es de rango completo.
La factorización
Esta factorización es usada posteriormente para resolver el
sistema
El esquema mostrado en la figura 1 puede ser usado para calcul=
ar
la inversa de la matriz A en problemas de tamańo moderado.
Una de las debilidades de la factorización LU es su
inestabilidad numérica; es por ello que en muchas aplicaciones es preferible
usar, a un mayor costo computacional, otro tipo de factorización, como la
factorización
Cuando la matriz A, de orden n, es simétrica y
definida positiva, una buena alternativa en términos de estabilidad, que
resulta menos costosa, es la factorización de Cholesky, la cual descompone =
a la
matriz A en el producto de una matriz triangular inferior L y=
su
traspuesta
Esta descomposición, en forma similar a la factorización LU=
,
utiliza el proceso basado en la eliminación Gaussiana para producir la matr=
iz
triangular L; la diferencia es que al hallar L ya se tiene qu=
e
Toda vez que la factorización de Cholesky posee característica=
s de
estabilidad, se reduce el costo computacional porque no se requiere pivoteo=
en
la eliminación Gaussiana. Así, en aplicaciones donde la matriz de los coefi=
cientes
es simétrica y definida positiva se prefiere la factorización de Cholesky p=
ara
resolver Ax=3Db. Más aún, en la práctica los algoritmos más eficient=
es no
usan la eliminación Gaussiana para hacer esta descomposición, sino que se
aprovechan de la estructura de la factorización, reduciendo aún más el costo
computacional.
La estructura de la factorización de Cholesky de una matriz
Una simple inspección de la figura 3 nos lleva al algoritmo de
Cholesky, el cual se basa en hallar los elementos de L en forma
iterativa. Conociendo el valor de
Cuando hablamos de matrices de
El método de los mínimos cuadrados es ampliamente usado en
aplicaciones de análisis de datos; en particular se propone en este trabajo=
, el
estudio de una de sus extensiones como lo es el método de los mínimos cuadr=
ados
reponderados, el cual es una extensión del conocido método de los mínimos
cuadrados ponderados que se describe en la próxima sección.
Esta extensión del método de los mínimos cuadrados, también
conocido como el método GSK, fue propuesto por Grizzie, Starmer y Koch en 1=
969,
para ser usado en la estimación de parámetros de un modelo lineal generaliz=
ado para
variables de respuesta con distribución multinomial, o de Poisson.
El método introduce una variante en el estudio de variables
categóricas, una vez que permite analizarlas a través de una función de
respuesta en lugar de sus frecuencias absolutas. Esta función de respuesta =
se
escogerá de acuerdo a la naturaleza de los datos y los objetivos que persig=
a la
investigación. Para mayor detalle sobre la evaluaión
y cálculo de estas funciones se recomienda Forthofer y Lehnen (1981).
De forma análoga al método de los mínimos cuadrados, este méto=
do
parte de un conjunto de observaciones
La escogencia de los pesos obedece a criterios para reducir o
amplificar la influencia de las observaciones en la recta de ajuste; en la
práctica se suele asociar a cada observación i un peso inversamente
proporcional a la variabilidad de yi, la cual se aproxima=
en
función de xi. En cualquier caso, si la suma de los pesos=
en
igual a 1, la solución del problema se basa en los cálculos que se presenta=
n en
la figura 4.
Ahora bien, tal como se ha dicho, las funciones=
de
peso pueden variar de acuerdo a una ponderación conveniente; de esta manera
pudiéramos tener funciones como:
1)&n=
bsp;
Ponderación uniforme:
2)&n=
bsp;
Ponderación Bisquare:
3)&n=
bsp;
Ponderación=
de
Huber:
4)&n=
bsp;
Ponderación
exponencial <=
!--[if gte msEquation 12]>
Existen otras funciones de ponderación que son =
más
apropiadas en ciertas aplicaciones, tales como la ponderación Gaussiana, el
spline cúbico, y otras que conllevan cálculos más pesados. En la práctica, =
las
más usadas son las ponderaciones de Bisquare y la de Huber. Para mayores
detalles sobre estas funciones refiérase a Torkzi (2011).
El método iterativo de los mínimos cuadrados re=
ponderados
(IRLS por sus siglas en inglés) resuelve en cada paso un problema usando
mínimos cuadrados ponderados, y resulta de mucho interés en aplicaciones do=
nde se
busca minimizar una norma p.
En particular el método es muy útil en el análi=
sis
de modelos lineales generalizados, en los cuales se estudia la relación:
El problema general que se resuelve es:
En este modelo,
Sin embargo, existen situaciones donde los erro=
res
se distribuyen por grupos y entre grupos difieren; en este caso el método de
los mínimos cuadrados ponderados es una mejor herramienta para estimar
También podemos encontrar situaciones que no se
relacionen con las
La figura 5 muestra los pasos básicos que conte=
mpla
el algoritmo iterativo de los mínimos cuadrados reponderados.
El valor inicial
En la práctica, todos estos métodos tienen una
forma matricial y la solución se halla resolviendo el correspondiente siste=
ma
de ecuaciones normales de <=
!--[if gte msEquation 12]>
En Ponsot (2011), el autor estudia la agregació=
n de
niveles de factores en un modelo logit binomial, y una de las propuestas es=
la
utilización del método de los mínimos cuadrados iterativos reponderados par=
a estimar
los parámetros a través de la función de verosimilitud. Las pruebas realiza=
das
con la implementación muestran el éxito de esa propuesta en la aplicación de
estudio.
Así como esta aplicación surgen otras en las cu=
ales
el método de los mínimos cuadrados iterativo ponderado resulta de utilidad.=
Es
por ello que proponemos en este trabajo, una implementación robusta de la
factorización de rango completo de Cholesky que mejore los tiempos de cómpu=
to
en las aproximaciones a los parámetros en cada iteración.
Una idea es implementar esquemas numéricos de la
factorización rectangular sobre la matriz de los coeficientes del sistema de
ecuaciones normales que se resuelve en cada iteración, en los cuales se evi=
ten
cálculos de inversas y la construcción de matrices de gran tamańo como es el
caso de
Una matriz A, deficiente en rango por
columnas, puede descomponerse en el producto de dos matrices rectangulares =
G
y U, tales que
ˇ
ˇ =
X: matriz de
diseńo de la experimentación, con valores no aleatorios.
ˇ
El vector de parámetros inicial
El algoritmo para los mínimos cuadrados iterati=
vo
reponderado utilizado fue implementado en Matlab y se presenta en la figura=
7.
Tal como puede observarse es una iteración que =
se
detiene cuando se satisface una tolerancia del error entre dos estimaciones
sucesivas. El código utiliza dos funciones auxiliares para aproximar la fun=
ción
Lo primero que notamos en el código presentado =
es
el cálculo de inversas generalizadas o pseudoinversa, en las líneas 3 y 12.=
Posteriormente
en la línea 7 se resuelve un sistema de ecuaciones con esta pseudoinversa. =
Se
propone en este trabajo sustituir estos cálculos por otros equivalentes
numéricamente que se ejecutan sobre la una factorización de rango completo =
de
las matrices involucradas. La resolución de sistemas de ecuaciones puede
hacerse en dos fases de menor costo con los factores de Cholesky.
Finalmente, una vez que en esta aplicación la
pseudoinversa de la matriz que se calcula en la línea 12 es requerida al fi=
nal
de este cálculo, se sugiere sustituir por una factorización rectangular de =
Cholesky
de rango completo, a fin de que solo sea necesario invertir el factor.
Los cambios sugeridos fueron implementados sobr=
e el
código original, usando la versión rectangular del algoritmo propuesto en la
figura 6 para hallar el factor de Cholesky de rango completo. La figura 8
muestra el código en Matlab con las modificaciones sugeridas.
Como puede observarse, en la versión propuesta =
se
sustituye el cálculo de las pseudoinversas en las líneas 3 y 12 del código
original. También se sustituye el producto con la pseudoinversa de la línea=
7
del código original, por la resolución de sistemas triangulares equivalente=
s con
el factor L y su traspuesta. Por último,
la pseudoinversa que se calcula en la línea 12 del código original, que es
necesaria para un futuro cálculo en esta aplicación, se ha sustituido por el
cálculo de la inversa del factor de Cholesky.
Se llevó a cabo una experimentación numérica
tomando matrices sintéticas X de n×m variando el número de fi=
las
entre n=3D500 y n=3D30.000, y manteniendo m =3D70. Se
seleccionó el vector y de respuestas aleatoriamente en cada
experimento. Los resultados compar=
ativos
entre la estimación hallada por el código original y la hallada por el códi=
go
modificado, así como también los tiempos de ejecución medidos se muestran e=
n la
figura 9.
Existen otras propuestas efectivas de carácter
general para acelerar el cálculo de estas factorizaciones, ver por ejemplo =
el
trabajo de Camero (2018). También pueden encontrarse eficientes
implementaciones en librerías desarrolladas en Matlab, Fortran, C, y C++, t=
ales
como las presentadas por Molt (2021) y Gustavson (2018). Sin embargo, para
aplicaciones específicas, como la que se estudia, la implementación propues=
ta
saca provecho de la estructura del problema y las características de la
solución buscada.
Tal como se puede observar en la figura 8, en
términos de precisión en la soluciones obtenidas, tanto el código original =
como
el modificado generan numéricamente la misma respuesta en la estimación de =
los
parámetros, en tanto que la versión modificada reduce significativamente el
tiempo de ejecución. Estos resultados son alentadores y se espera que puedan
corroborarse con datos experimentales reales.
ˇ =
Hemos
presentado en este trabajo una propuesta algorítmica para la resolución de =
los
mínimos cuadrados iterativos reponderados, usando una factorización de rango
completo de Cholesky cuando la matriz es rectangular o deficiente en rango.=
Se
resaltan las ventajas de su aplicación en algoritmos estadísticos que requi=
eren
la resolución de sistemas de ecuaciones normales donde la matriz de los
coeficientes es rectangular o posiblemente con columnas dependientes; situa=
ción
que se presenta con frecuencia en análisis de datos categóricos.
ˇ =
En
particular, como un estudio de caso, se analizó la factorización de Cholesk=
y de
rango completo en una aplicación de los mínimos cuadrados iterativos
reponderados para la estimación estadística máximo verosímil de parámetros =
en
un modelo logit binomial. Los resultados sugieren que puede mantenerse la
precisión en la solución obtenida al mismo tiempo que se reduce
considerablemente el tiempo de cómputo.
ˇ =
Tomando
en cuenta que el algoritmo propuesto está basado mayormente en operaciones =
con
matrices y vectores, los cuales son altamente paralelizables, el próximo pa=
so
en esta dirección será la evaluación numérico-computacional de la
implementación de esta descomposición en máquinas de alto rendimiento o con
múltiples unidades de procesamiento.
<=
span
lang=3DEN-US style=3D'font-size:12.0pt;line-height:115%;font-family:"Times =
New Roman",serif;
color:black;mso-ansi-language:EN-US'>Referencias bibliográficas
Bjorck, A. (1996).
Numerical Method for Least Squares Problems, Society for Indus=
trial
and Applied Mathematics; 1a. Edición.
Camero, C. (2018). Simp=
le,
Fast and Practicable Algorithms for Cholesky, LU and QR Decomposition using
fast rectangular matrix multiplications, Department of Computer Science and
Electronics, Universidad de Cantabria, UNICAN, Spain.
Cantó, R., Ricarte B.,
Urbano. A.M. (2009).
Full rank factorization and Flanders Theorem, Electronic Journal of Linear
Algebra, vol. 18, 352-363.
Cantó, R., Pel=
áez
M., Urbano. A.=
M.
(2015). Full rank Cholesky factorization for rank deficient matrices. Appli=
ed
Mathematics Letters, 40, 17-22.
Chen, J. (1993). Iterative Weighted Le=
ast
Square Estimators, The
Annals of Statistics, 21(2).
Davidian, M. a=
nd Carroll,
R. J. (1987). Variance function estimation. J. Amer. Statist. Assoc.=
Vol.
82, No. 400,
pp. 1079-1091.
Datta, B.N (20=
10).
Numerical Linear Algebra and Applications. Society for Industr=
ial
and Applied Mathematics; 2a. Edición.
Forthofer R.N.,
Lehnen R.G. (1981). One Response and Two Factor Variables. In: Public Progr=
am
Analysis. Springer, Boston, MA. https://doi.org/10.1007/978-1-4684-6683-6_4.
Fuller, W. A. =
and
Rao, J. N. K. (1978). Estimation for a linear regression model with unknown
diagonal covariance matrix. Ann. Statist., 6, 1149 - 1158.
Gustavson, F.,
Wasniewski, J., Dongarra, J., Langou, J. (2009). Rectangular Dull Packed Fo=
rmat
for Choleskys Algorithm: Factorization, Solution and Inversion ACM Transac=
tion
on Mathematical Software.
Molt, J. (2021=
).
Cholesky decomposition. Recuperado August 6, 2021 de MATLAB Central File Ex=
change.
htt=
ps://www.mathworks.com/matlabcentral/fileexchange/71304-cholesky-decomposit=
ion.
Ponsot, E.
(2011). Estudio de la agregación de niveles de=
los
factores en el modelo logit binomial. Tesis Doctoral. Instituto de Estadíst=
ica
Aplicada, Universidad de los Andes, Venezuela.
Torkzi, K. (20=
11)
Aproximaciones por mínimos cuadrados ponderados en el entorno de
multirresolución por valores puntuales. Departamento de Matemática Aplicada=
y
Estadística, Universidad Politécnica de Cartagena, Espańa.
Trefethen, L. N. & Bau
III, D. (1997). Numerical linear algebra. Siam, Philadelphia.
PARA
CITAR EL ARTÍCULO INDEXADO.
Castillo Marrero, Z. N., Ponsot
Balaguer, E. A., José Camacho, F., & León Sánchez, M. V.
. (2021). Usando la Factorización de Rango Completo de Cholesky en la Solución de los Mínimos Cuadrados Ponderados . AlfaPublicaciones,
3(3.1). https://doi.org/1=
0.33262/ap.v3i3.1.79
El
artículo que se publica es de exclusiva responsabilidad de los autores y no
necesariamente reflejan el pensamiento de la Revista Alfa Publicaciones.
El artículo qu=
eda
en propiedad de la revista y, por tanto, su publicación parcial y/o total en
otro medio tiene que ser autorizado por el director de la Revista Alfa Publicaciones.
=
=
[1] Escue=
la
Superior Politécnica de Chimborazo, Facultad de Ciencias, Escuela de Matemá=
tica,
Grupo CIDED, Riobamba, Ecuador, zenaida.castillo@espoch.edu.ec,
https://orcid.org/0000-0002-4424-8652.
[2] <= span style=3D'font-family:"Times New Roman",serif'>Data Sci= ence Consulting/Universidad de los Andes, Escuela de Estadística. Mérida, Venezuela, ernesto.pb@gmail.com, https://or= cid.org/0000-0001-5221-1799.
[3] Universidad de Investigación=
de
Tecnología Experimental Yachay, Escuela
de Ciencias Matemáticas y Computacionales, Urcuquí,060150, Ecuador, =
fcamacho@yachaytech.edu.ec,
https://orcid.org/=
0000-0001=
-7802-5687
[4] Escue=
la
Superior Politécnica de Chimborazo, Facultad de Ciencias, Escuela de Matemá=
tica,
Riobamba, Ecuador, maria.leon@espoch.edu.ec, https://orcid.org/0000-0003-06=
12-7478.
Alpha publicaciones<=
o:p>