Aproximación de bajo rango

En matemáticas, la aproximación de bajo rango se refiere al proceso de aproximar una matriz dada por una matriz de rango inferior. Más precisamente, es un problema de minimización , en el que la función de costo mide el ajuste entre una matriz dada (los datos) y una matriz de aproximación (la variable de optimización), sujeta a una restricción de que la matriz de aproximación tiene un rango reducido . El problema se utiliza para el modelado matemático y la compresión de datos . La restricción de rango está relacionada con una restricción en la complejidad de un modelo que se ajusta a los datos. En las aplicaciones, a menudo hay otras restricciones en la matriz de aproximación aparte de la restricción de rango, por ejemplo, la no negatividad y la estructura de Hankel .

La aproximación de bajo rango está estrechamente relacionada con numerosas otras técnicas, incluido el análisis de componentes principales , el análisis factorial , los mínimos cuadrados totales , el análisis semántico latente , la regresión ortogonal y la descomposición en modos dinámicos .

Dado

• especificación de estructura ,${\displaystyle {\mathcal {S}}:\mathbb {R} ^{n_{p}}\to \mathbb {R} ^{m\times n}}$
• vector de parámetros de estructura ,${\displaystyle p\in \mathbb {R} ^{n_{p}}}$
• norma , y${\estilo de visualización \|\cdot \|}$
• rango deseado ,${\estilo de visualización r}$

${\displaystyle {\text{minimizar}}\quad {\text{sobre }}{\widehat {p}}\quad \|p-{\widehat {p}}\|\quad {\text{sujeto a}}\quad \operatorname {rango} {\big (}{\mathcal {S}}({\widehat {p}}){\big )}\leq r.}$

• Identificación de sistema lineal , en cuyo caso la matriz de aproximación está estructurada de Hankel .
• Aprendizaje automático , en cuyo caso la matriz de aproximación está estructurada de forma no lineal.
• Sistemas de recomendación , en cuyos casos la matriz de datos tiene valores faltantes y la aproximación es categórica .
• Completitud de la matriz de distancia , en cuyo caso hay una restricción de definición positiva.
• Procesamiento del lenguaje natural , en cuyo caso la aproximación es no negativa .
• Álgebra computacional , en cuyo caso la aproximación es estructurada de Sylvester .

El problema no estructurado con ajuste medido por la norma de Frobenius , es decir,

${\displaystyle {\text{minimizar}}\quad {\text{sobre }}{\widehat {D}}\quad \|D-{\widehat {D}}\|_{\text{F}}\quad {\text{sujeto a}}\quad \operatorname {rango} {\big (}{\widehat {D}}{\big )}\leq r}$

tiene una solución analítica en términos de la descomposición en valores singulares de la matriz de datos. El resultado se conoce como lema de aproximación matricial o teorema de Eckart-Young-Mirsky . Este problema fue resuelto originalmente por Erhard Schmidt ^[1] en el contexto de dimensión infinita de operadores integrales (aunque sus métodos se generalizan fácilmente a operadores compactos arbitrarios en espacios de Hilbert) y luego redescubierto por C. Eckart y G. Young . ^[2] L. Mirsky generalizó el resultado a normas unitariamente invariantes arbitrarias. ^[3] Sea

${\displaystyle D=U\Sigma V^{\top }\in \mathbb {R} ^{m\times n},\quad m\geq n}$

sea ​​la descomposición en valores singulares de , donde es la matriz diagonal rectangular con los valores singulares . Para un determinado , particione , , y de la siguiente manera: ${\estilo de visualización D}$${\displaystyle \Sigma =:\nombre del operador {diag} (\sigma _{1},\ldots ,\sigma _{m})}$${\estilo de visualización m\veces m}$${\displaystyle \sigma _{1}\geq \ldots \geq \sigma _{m}}$${\displaystyle r\en \{1,\puntos ,m-1\}}$${\estilo de visualización U}$${\estilo de visualización \Sigma}$${\estilo de visualización V}$

${\displaystyle U=:{\begin{bmatrix}U_{1}&U_{2}\end{bmatrix}},\quad \Sigma =:{\begin{bmatrix}\Sigma _{1}&0\\0&\Sigma _{2}\end{bmatrix}},\quad {\text{y}}\quad V=:{\begin{bmatrix}V_{1}&V_{2}\end{bmatrix}},}$

donde es , es , y es . Entonces la matriz de rango, obtenida a partir de la descomposición en valores singulares truncados $Estilo de visualización U_{1}$${\displaystyle m\times r}$${\displaystyle \Sigma _{1}}$${\displaystyle r\times r}$${\displaystyle V_{1}}$${\displaystyle r\times n}$${\displaystyle r}$

${\displaystyle {\widehat {D}}^{*}=U_{1}\Sigma _{1}V_{1}^{\top },}$

es tal que

${\displaystyle \|D-{\widehat {D}}^{*}\|_{\text{F}}=\min _{\operatorname {rank} ({\widehat {D}})\leq r}\|D-{\widehat {D}}\|_{\text{F}}={\sqrt {\sigma _{r+1}^{2}+\cdots +\sigma _{m}^{2}}}.}$

El minimizador es único si y sólo si . ${\displaystyle {\widehat {D}}^{*}}$${\displaystyle \sigma _{r+1}\neq \sigma _{r}}$

Sea una matriz real (posiblemente rectangular) con . Supóngase que ${\displaystyle A\in \mathbb {R} ^{m\times n}}$${\displaystyle m\leq n}$

${\displaystyle A=U\Sigma V^{\top }}$

es la descomposición en valores singulares de . Recuerde que y son matrices ortogonales , y es una matriz diagonal con entradas tales que . ${\displaystyle A}$${\displaystyle U}$${\displaystyle V}$${\displaystyle \Sigma }$${\displaystyle m\times n}$ ${\displaystyle (\sigma _{1},\sigma _{2},\cdots ,\sigma _{m})}$${\displaystyle \sigma _{1}\geq \sigma _{2}\geq \cdots \geq \sigma _{m}\geq 0}$

Afirmamos que la mejor aproximación de rango a en la norma espectral, denotada por , viene dada por ${\displaystyle k}$${\displaystyle A}$${\displaystyle \|\cdot \|_{2}}$

${\displaystyle A_{k}:=\sum _{i=1}^{k}\sigma _{i}u_{i}v_{i}^{\top }}$

donde y denotan la ésima columna de y , respectivamente. ${\displaystyle u_{i}}$${\displaystyle v_{i}}$${\displaystyle i}$${\displaystyle U}$${\displaystyle V}$

En primer lugar, tenga en cuenta que tenemos

${\displaystyle \|A-A_{k}\|_{2}=\left\|\sum _{i=1}^{\color {red}{n}}\sigma _{i}u_{i}v_{i}^{\top }-\sum _{i=1}^{\color {red}{k}}\sigma _{i}u_{i}v_{i}^{\top }\right\|_{2}=\left\|\sum _{i=\color {red}{k+1}}^{n}\sigma _{i}u_{i}v_{i}^{\top }\right\|_{2}=\sigma _{k+1}}$

Por lo tanto, necesitamos demostrar que si donde y tienen columnas entonces . ${\displaystyle B_{k}=XY^{\top }}$${\displaystyle X}$${\displaystyle Y}$${\displaystyle k}$${\displaystyle \|A-A_{k}\|_{2}=\sigma _{k+1}\leq \|A-B_{k}\|_{2}}$

Dado que tiene columnas, entonces debe haber una combinación lineal no trivial de las primeras columnas de , es decir, ${\displaystyle Y}$${\displaystyle k}$${\displaystyle k+1}$${\displaystyle V}$

${\displaystyle w=\gamma _{1}v_{1}+\cdots +\gamma _{k+1}v_{k+1},}$

de modo que . Sin pérdida de generalidad, podemos escalar de modo que o (equivalentemente) . Por lo tanto, ${\displaystyle Y^{\top }w=0}$${\displaystyle w}$${\displaystyle \|w\|_{2}=1}$${\displaystyle \gamma _{1}^{2}+\cdots +\gamma _{k+1}^{2}=1}$

${\displaystyle \|A-B_{k}\|_{2}^{2}\geq \|(A-B_{k})w\|_{2}^{2}=\|Aw\|_{2}^{2}=\gamma _{1}^{2}\sigma _{1}^{2}+\cdots +\gamma _{k+1}^{2}\sigma _{k+1}^{2}\geq \sigma _{k+1}^{2}.}$

El resultado se obtiene tomando la raíz cuadrada de ambos lados de la desigualdad anterior.

Sea una matriz real (posiblemente rectangular) con . Supóngase que ${\displaystyle A\in \mathbb {R} ^{m\times n}}$${\displaystyle m\leq n}$

${\displaystyle A=U\Sigma V^{\top }}$

es la descomposición en valor singular de . ${\displaystyle A}$

Afirmamos que la mejor aproximación de rango a en la norma de Frobenius, denotada por , viene dada por ${\displaystyle k}$${\displaystyle A}$${\displaystyle \|\cdot \|_{F}}$

${\displaystyle A_{k}=\sum _{i=1}^{k}\sigma _{i}u_{i}v_{i}^{\top }}$

donde y denotan la ésima columna de y , respectivamente. ${\displaystyle u_{i}}$${\displaystyle v_{i}}$${\displaystyle i}$${\displaystyle U}$${\displaystyle V}$

En primer lugar, tenga en cuenta que tenemos

${\displaystyle \|A-A_{k}\|_{F}^{2}=\left\|\sum _{i=k+1}^{n}\sigma _{i}u_{i}v_{i}^{\top }\right\|_{F}^{2}=\sum _{i=k+1}^{n}\sigma _{i}^{2}}$

Por lo tanto, necesitamos demostrar que si donde y tienen columnas entonces ${\displaystyle B_{k}=XY^{\top }}$${\displaystyle X}$${\displaystyle Y}$${\displaystyle k}$

${\displaystyle \|A-A_{k}\|_{F}^{2}=\sum _{i=k+1}^{n}\sigma _{i}^{2}\leq \|A-B_{k}\|_{F}^{2}.}$

Por la desigualdad triangular con la norma espectral, si entonces . Supóngase que y denotan respectivamente la aproximación de rango a y por el método SVD descrito anteriormente. Entonces, para cualquier${\displaystyle A=A'+A''}$${\displaystyle \sigma _{1}(A)\leq \sigma _{1}(A')+\sigma _{1}(A'')}$${\displaystyle A'_{k}}$${\displaystyle A''_{k}}$${\displaystyle k}$${\displaystyle A'}$${\displaystyle A''}$${\displaystyle i,j\geq 1}$

${\displaystyle {\begin{aligned}\sigma _{i}(A')+\sigma _{j}(A'')&=\sigma _{1}(A'-A'_{i-1})+\sigma _{1}(A''-A''_{j-1})\\&\geq \sigma _{1}(A-A'_{i-1}-A''_{j-1})\\&\geq \sigma _{1}(A-A_{i+j-2})\qquad ({\text{since }}{\rm {rank}}(A'_{i-1}+A''_{j-1})\leq i+j-2))\\&=\sigma _{i+j-1}(A).\end{aligned}}}$

Desde , cuando y concluimos que para${\displaystyle \sigma _{k+1}(B_{k})=0}$${\displaystyle A'=A-B_{k}}$${\displaystyle A''=B_{k}}$${\displaystyle i\geq 1,j=k+1}$

${\displaystyle \sigma _{i}(A-B_{k})\geq \sigma _{k+i}(A).}$

Por lo tanto,

${\displaystyle \|A-B_{k}\|_{F}^{2}=\sum _{i=1}^{n}\sigma _{i}(A-B_{k})^{2}\geq \sum _{i=k+1}^{n}\sigma _{i}(A)^{2}=\|A-A_{k}\|_{F}^{2},}$

según sea necesario.

La norma de Frobenius pondera de manera uniforme todos los elementos del error de aproximación . El conocimiento previo sobre la distribución de los errores se puede tener en cuenta considerando el problema de aproximación ponderada de bajo rango. ${\displaystyle D-{\widehat {D}}}$

${\displaystyle {\text{minimize}}\quad {\text{over }}{\widehat {D}}\quad \operatorname {vec} (D-{\widehat {D}})^{\top }W\operatorname {vec} (D-{\widehat {D}})\quad {\text{subject to}}\quad \operatorname {rank} ({\widehat {D}})\leq r,}$

donde vectoriza la matriz columna por columna y es una matriz de peso (semi)definido positivo dado. ${\displaystyle {\text{vec}}(A)}$ ${\displaystyle A}$${\displaystyle W}$

El problema general de aproximación ponderada de bajo rango no admite una solución analítica en términos de la descomposición en valores singulares y se resuelve mediante métodos de optimización local, que no ofrecen ninguna garantía de que se encuentre una solución globalmente óptima.

En el caso de pesos no correlacionados, el problema de aproximación ponderada de rango bajo también se puede formular de esta manera: ^{[4] [5]} para una matriz no negativa y una matriz que queremos minimizar sobre matrices, , de rango como máximo . ${\displaystyle W}$${\displaystyle A}$${\displaystyle \sum _{i,j}(W_{i,j}(A_{i,j}-B_{i,j}))^{2}}$${\displaystyle B}$${\displaystyle r}$

Sea . Para , el algoritmo más rápido se ejecuta en el tiempo. ^{[6] [7]} Una de las ideas importantes que se han utilizado se llama Oblivious Subspace Embedding (OSE), propuesta por primera vez por Sarlos. ^[8]${\displaystyle \|A\|_{p}=\left(\sum _{i,j}|A_{i,j}^{p}|\right)^{1/p}}$${\displaystyle p=2}$${\displaystyle nnz(A)+n\cdot poly(k/\epsilon )}$

Para , se sabe que esta norma L1 de entrada por entrada es más robusta que la norma de Frobenius en presencia de valores atípicos y está indicada en modelos donde los supuestos gaussianos sobre el ruido pueden no aplicarse. Es natural buscar minimizar . ^[9] Para y , hay algunos algoritmos con garantías demostrables. ^{[10] [11]}${\displaystyle p=1}$${\displaystyle \|B-A\|_{1}}$${\displaystyle p=0}$${\displaystyle p\geq 1}$

Sean y dos conjuntos de puntos en un espacio métrico arbitrario. Sea la matriz donde . Estas matrices de distancias se calculan comúnmente en paquetes de software y tienen aplicaciones para el aprendizaje de variedades de imágenes, reconocimiento de escritura a mano y desdoblamiento multidimensional. En un intento por reducir su tamaño de descripción, ^{[12] [13]} se puede estudiar la aproximación de bajo rango de dichas matrices. ${\displaystyle P=\{p_{1},\ldots ,p_{m}\}}$${\displaystyle Q=\{q_{1},\ldots ,q_{n}\}}$${\displaystyle A}$${\displaystyle m\times n}$${\displaystyle A_{i,j}=dist(p_{i},q_{i})}$

Los problemas de aproximación de bajo rango en el entorno distribuido y de transmisión se han considerado en ^{[14] .}

Usando las equivalencias

${\displaystyle \operatorname {rank} ({\widehat {D}})\leq r\quad \iff \quad {\text{there are }}P\in \mathbb {R} ^{m\times r}{\text{ and }}L\in \mathbb {R} ^{r\times n}{\text{ such that }}{\widehat {D}}=PL}$

y

${\displaystyle \operatorname {rank} ({\widehat {D}})\leq r\quad \iff \quad {\text{there is full row rank }}R\in \mathbb {R} ^{m-r\times m}{\text{ such that }}R{\widehat {D}}=0}$

El problema de aproximación de rango bajo ponderado se vuelve equivalente a los problemas de optimización de parámetros.

${\displaystyle {\text{minimize}}\quad {\text{over }}{\widehat {D}},P{\text{ and }}L\quad \operatorname {vec} ^{\top }(D-{\widehat {D}})W\operatorname {vec} (D-{\widehat {D}})\quad {\text{subject to}}\quad {\widehat {D}}=PL}$

y

${\displaystyle {\text{minimize}}\quad {\text{over }}{\widehat {D}}{\text{ and }}R\quad \operatorname {vec} ^{\top }(D-{\widehat {D}})W\operatorname {vec} (D-{\widehat {D}})\quad {\text{subject to}}\quad R{\widehat {D}}=0\quad {\text{and}}\quad RR^{\top }=I_{r},}$

donde es la matriz identidad de tamaño . ${\displaystyle I_{r}}$${\displaystyle r}$

La representación en imagen de la restricción de rango sugiere un método de optimización de parámetros en el que la función de costo se minimiza alternativamente sobre una de las variables ( o ) con la otra fija. Aunque la minimización simultánea sobre ambas y es un problema de optimización biconvexa difícil , la minimización sobre una de las variables por sí sola es un problema de mínimos cuadrados lineales y se puede resolver de manera global y eficiente. ${\displaystyle P}$${\displaystyle L}$${\displaystyle P}$${\displaystyle L}$

El algoritmo de optimización resultante (llamado proyecciones alternas) es convergente globalmente con una tasa de convergencia lineal hacia una solución localmente óptima del problema de aproximación de rango bajo ponderado. Se debe indicar el valor inicial del parámetro (o ). La iteración se detiene cuando se satisface una condición de convergencia definida por el usuario. ${\displaystyle P}$${\displaystyle L}$

Implementación de Matlab del algoritmo de proyecciones alternas para aproximación ponderada de rango bajo:

función [dh, f] = wlra_ap ( d, w, p, tol, maxiter ) [ m , n ] = tamaño ( d ); r = tamaño ( p , 2 ); f = inf ; para i = 2 : maxiter % minimización sobre L bp = kron ( ojo ( n ), p ); vl = ( bp ' * w * bp ) \ bp ' * w * d (:); l = remodelar ( vl , r , n ); % minimización sobre P bl = kron ( l ' , ojo ( m )); vp = ( bl ' * w * bl ) \ bl ' * w * d (:); p = remodelar ( vp , m , r ); % comprobar condición de salida dh = p * l ; dd = d - dh ; f ( i ) = dd (:) ' * w * dd (:); si abs ( f ( i - 1 ) - f ( i )) < tol , break , fin finpara 
          
   
    
       
                
         
    
       
                
         
    
             
          
             

El algoritmo de proyecciones alternas aprovecha el hecho de que el problema de aproximación de rango bajo, parametrizado en forma de imagen, es bilineal en las variables o . La naturaleza bilineal del problema se utiliza de manera efectiva en un enfoque alternativo, llamado proyecciones variables. ^[15]${\displaystyle P}$${\displaystyle L}$

Consideremos nuevamente el problema de aproximación de rango bajo ponderado, parametrizado en forma de imagen. La minimización con respecto a la variable (un problema de mínimos cuadrados lineal) conduce a la expresión en forma cerrada del error de aproximación como una función de${\displaystyle L}$${\displaystyle P}$

${\displaystyle f(P)={\sqrt {\operatorname {vec} ^{\top }(D){\Big (}W-W(I_{n}\otimes P){\big (}(I_{n}\otimes P)^{\top }W(I_{n}\otimes P){\big )}^{-1}(I_{n}\otimes P)^{\top }W{\Big )}\operatorname {vec} (D)}}.}$

Por lo tanto, el problema original es equivalente al problema de mínimos cuadrados no lineales de minimizar con respecto a . Para este propósito, se pueden utilizar métodos de optimización estándar, por ejemplo, el algoritmo de Levenberg-Marquardt .${\displaystyle f(P)}$${\displaystyle P}$

Implementación de Matlab del algoritmo de proyecciones de variables para la aproximación ponderada de rango bajo:

función [dh, f] = wlra_varpro ( d, w, p, tol, maxiter ) prob = optimset (); prob . solver = 'lsqnonlin' ; prob . opciones = optimset ( 'MaxIter' , maxiter , 'TolFun' , tol ); prob . x0 = p ; prob . objetivo = @ ( p ) costo_diversión ( p , d , w ); [ p , f ] = lsqnonlin ( prob ); [ f , vl ] = costo_diversión ( p , d , w ); dh = p * remodelar ( vl , tamaño ( p , 2 ), tamaño ( d , 2 )); 
     
      
        
     
      
        

función [f, vl] = costo_diversión ( p, d, w ) bp = kron ( ojo ( tamaño ( d , 2 )), p ); vl = ( bp ' * w * bp ) \ bp ' * w * d (:); f = d (:) ' * w * ( d (:) - bp * vl ); 
    
            
          

El enfoque de proyecciones de variables también se puede aplicar a problemas de aproximación de bajo rango parametrizados en la forma kernel. El método es efectivo cuando el número de variables eliminadas es mucho mayor que el número de variables de optimización que quedan en la etapa de minimización de mínimos cuadrados no lineales. Tales problemas ocurren en la identificación de sistemas, parametrizados en la forma kernel, donde las variables eliminadas son la trayectoria de aproximación y las variables restantes son los parámetros del modelo. En el contexto de sistemas lineales invariantes en el tiempo , el paso de eliminación es equivalente al suavizado de Kalman .

Por lo general, queremos que nuestra nueva solución no solo sea de bajo rango, sino que también satisfaga otras restricciones convexas debido a los requisitos de la aplicación. El problema que nos interesa sería el siguiente:

${\displaystyle {\text{minimize}}\quad {\text{over }}{\widehat {p}}\quad \|p-{\widehat {p}}\|\quad {\text{subject to}}\quad \operatorname {rank} {\big (}{\mathcal {S}}({\widehat {p}}){\big )}\leq r{\text{ and }}g({\widehat {p}})\leq 0}$

Este problema tiene muchas aplicaciones en el mundo real, incluida la recuperación de una buena solución a partir de una relajación inexacta (programación semidefinida). Si la restricción adicional es lineal, como si requiriéramos que todos los elementos sean no negativos, el problema se denomina aproximación de bajo rango estructurada. ^[16] La forma más general se denomina aproximación de bajo rango con restricción convexa. ${\displaystyle g({\widehat {p}})\leq 0}$

Este problema es útil para resolver muchos problemas. Sin embargo, es un desafío debido a la combinación de las restricciones convexas y no convexas (de bajo rango). Se desarrollaron diferentes técnicas basadas en diferentes realizaciones de . Sin embargo, el método de dirección alternada de multiplicadores (ADMM) se puede aplicar para resolver el problema no convexo con función objetivo convexa, restricciones de rango y otras restricciones convexas, ^[17] y, por lo tanto, es adecuado para resolver nuestro problema anterior. Además, a diferencia de los problemas no convexos generales, ADMM garantizará la convergencia de una solución factible siempre que su variable dual converja en las iteraciones. ${\displaystyle g({\widehat {p}})\leq 0}$

• La aproximación de la matriz CUR se realiza a partir de las filas y columnas de la matriz original.

• MT Chu, RE Funderlic, RJ Plemmons, Aproximación estructurada de bajo rango, Álgebra lineal y sus aplicaciones, Volumen 366, 1 de junio de 2003, páginas 157–172 doi :10.1016/S0024-3795(02)00505-0

• Paquete C++ para aproximación estructurada de rango bajo