{"id":301,"date":"2020-07-13T19:12:01","date_gmt":"2020-07-13T17:12:01","guid":{"rendered":"https:\/\/intsight.com\/?p=301"},"modified":"2021-01-21T17:00:48","modified_gmt":"2021-01-21T16:00:48","slug":"cholesky","status":"publish","type":"post","link":"https:\/\/intsight.com\/index.php\/2020\/07\/13\/cholesky\/","title":{"rendered":"Cholesky"},"content":{"rendered":"

Esta va a ser una entrada<\/span> \u00abutilitaria\u00bb: poca explicaci\u00f3n, pero aportando c\u00f3digo fuente, por si alguien lo necesita. \u00bfRecuerda la entrada en la que present\u00e9 la descomposici\u00f3n de Cholesky<\/a>? En aquel momento, no inclu\u00ed el algoritmo de la descomposici\u00f3n, porque quer\u00eda experimentar un poco con la implementaci\u00f3n. Ya lo he hecho, y estoy m\u00e1s o menos satisfecho con el resultado.<\/p>\n

Antes de ver el c\u00f3digo, le explico el contexto. Este es un m\u00e9todo de una estructura Matrix<\/code>, que encapsula solamente una matriz bidimensional de valores flotantes (campo values<\/code>). Dentro de la estructura se definen todos los operadores y m\u00e9todos que pod\u00e9is imaginar. En paralelo, defino otra estructura, LowerMatrix<\/code>, para representar matrices triangulares inferiores. No hay relaciones de herencia, al tratarse de estructura, pero la clase de matrices triangulares permite definir m\u00e9todos y operadores m\u00e1s eficientes. Es importante tener en cuenta que LowerMatrix<\/code> gasta exactamente la misma memoria que Matrix<\/code>. La ventaja est\u00e1 en esos m\u00e9todos que evitan procesar las dos mitades de la matriz.<\/p>\n

Tambi\u00e9n hace falta saber que he definido operadores de conversi\u00f3n impl\u00edcitos que transforman un array bidimensional en uno u otro tipo de matrices. Este operador es el que me permite evitar un constructor expl\u00edcito para devolver un valor en el m\u00e9todo:<\/p>\n

public unsafe LowerMatrix Cholesky()\n{\n    int n = Rows;\n    double[,] dest = new double[n, n];\n    double[,] src = values;\n    double* tmp = stackalloc double[n + n];\n\n    \/\/ First column is special.\n    double ajj = src[0, 0];\n    if (ajj <= 0)\n    {\n        dest[0, 0] = double.NaN;\n        return dest;\n    }\n    dest[0, 0] = ajj = Math.Sqrt(ajj);\n    double r = 1 \/ ajj;\n    n--;\n    for (int i = 1; i <= n; i++)\n        dest[i, 0] = src[i, 0] * r;\n    for (int j = 1; j <= n; j++)\n    {\n        \/\/ Compute the diagonal cell.\n        double v = 0.0;\n        for (int i = 0; i < j; i++)\n        {\n            double a = dest[j, i];\n            v += a * a;\n        }\n        ajj = src[j, j] - v;\n        if (ajj <= 0)\n        {\n            dest[j, j] = double.NaN;\n            return dest;\n        }\n        dest[j, j] = ajj = Math.Sqrt(ajj);\n\n        \/\/ Compute the other cells of column J.\n        if (j < n)\n        {\n            r = 1 \/ ajj;\n            for (int i = 0; i < j; i++)\n                tmp[i] = dest[j, i];\n            for (int i = j; i < n; i++)\n            {\n                v = 0.0;\n                for (int k = 0; k < j; k++)\n                    v += dest[i + 1, k] * tmp[k];\n                tmp[i] = v;\n            }\n            for (int i = j; i < n; i++)\n                dest[i + 1, j] = (src[i + 1, j] - tmp[i]) * r;\n        }\n    }\n    return dest;\n}<\/pre>\n
Sobre el algoritmo, en s\u00ed: el algoritmo de Cholesky puede fallar cuando la matriz de origen no es positiva semidefinida. Mi implementaci\u00f3n detecta ese caso al calcular las ra\u00edces cuadradas… y se limita a parar, poniendo un NaN<\/code> en la celda diagonal donde se ha detectado el problema. Esto significa que el m\u00e9todo, en la pr\u00e1ctica, asume que la matriz es positiva semidefinida. En mi c\u00f3digo, tengo un segundo m\u00e9todo, TryCholesky<\/code>, que devuelve un valor l\u00f3gico para ver si la conversi\u00f3n fue posible, y retorna la matriz transformada como par\u00e1metro de salida.<\/p>\n
Desde el punto de vista de un programador de C#, el \u00fanico detalle interesante es el uso de stackalloc<\/strong> para reservar un array de memoria en la pila, en vez de usar memoria din\u00e1mica. Esto es lo que obliga a declarar el m\u00e9todo con unsafe<\/strong>.<\/p>\n
En rendimiento, el m\u00e9todo es m\u00e1s r\u00e1pido que la \u00abversi\u00f3n base\u00bb de la librer\u00eda que he visto que es m\u00e1s r\u00e1pida usando s\u00f3lo C# (usando Intel MKL, no hay competencia posible). Me refiero a la versi\u00f3n base porque, para matrices grandes, las librer\u00edas serias suelen dividir la matriz en bloques que se pueden procesar en paralelo. Este c\u00f3digo, como ve, no usa threads<\/em>, instrucciones SIMD y s\u00f3lo utiliza punteros para la cach\u00e9. En menos palabras: todo es mejorable.<\/p>\n","protected":false},"excerpt":{"rendered":"
Esta va a ser una entrada \u00abutilitaria\u00bb: poca explicaci\u00f3n, pero aportando c\u00f3digo fuente, por si alguien lo necesita. \u00bfRecuerda la entrada en la que present\u00e9 la descomposici\u00f3n de Cholesky? En aquel momento, no inclu\u00ed el algoritmo de la descomposici\u00f3n, porque quer\u00eda experimentar un poco con la implementaci\u00f3n. Ya lo he hecho, y estoy m\u00e1s o […]<\/p>\n","protected":false},"author":1,"featured_media":302,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"jetpack_post_was_ever_published":false,"_jetpack_newsletter_access":"","_jetpack_dont_email_post_to_subs":false,"_jetpack_newsletter_tier_id":0,"_jetpack_memberships_contains_paywalled_content":false,"_jetpack_memberships_contains_paid_content":false,"footnotes":""},"categories":[4,7],"tags":[15,17,19],"class_list":["post-301","post","type-post","status-publish","format-standard","has-post-thumbnail","hentry","category-c","category-fintech","tag-algorithms","tag-cholesky","tag-linear-algebra"],"jetpack_featured_media_url":"https:\/\/i0.wp.com\/intsight.com\/wp-content\/uploads\/2020\/07\/block.png?fit=350%2C350&ssl=1","jetpack_sharing_enabled":true,"_links":{"self":[{"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/posts\/301","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/comments?post=301"}],"version-history":[{"count":8,"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/posts\/301\/revisions"}],"predecessor-version":[{"id":558,"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/posts\/301\/revisions\/558"}],"wp:featuredmedia":[{"embeddable":true,"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/media\/302"}],"wp:attachment":[{"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/media?parent=301"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/categories?post=301"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/tags?post=301"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}