C#，数值计算，矩阵的乔莱斯基分解（Cholesky decomposition）算法与源代码

一、安德烈·路易斯·乔尔斯基

安德烈·路易斯·乔尔斯基出生于法国波尔多以北的查伦特斯海域的蒙古扬。他在波尔多参加了Lycée e，并于1892年11月14日获得学士学位的第一部分，于1893年7月24日获得第二部分。1895年10月15日，乔尔斯基进入莱科尔理工学院，在当年223名入学学生中排名第88位。他在莱科尔理工学院的教授包括卡米尔·乔丹和发现放射性的著名物理学家亨利·贝克勒尔。在成功的两年后，他于1897年参加了莱科尔理工学院的期末考试。在222名学生中，他提高了自己的地位，在这些考试中排名第38位。随后，他加入军队，成为少尉，并从1897年10月开始在炮兵学校学习。他在1899年完成了学业，并保持了稳步的进步，因为现在他在那一年获得资格的86名学生中排名第五。

突尼斯当时仍然是奥斯曼帝国的正式一部分，但在1883年突尼斯和法国政府签署的《马尔萨公约》之后，它一直是法国的保护国。法国人稳定了经济，建立了现代通讯。从1902年1月到6月，乔尔斯基在突尼斯执行了一项任务，然后从1902年11月到1903年5月，他被派去执行第二项任务。1903年12月31日，乔尔斯基开始在阿尔及利亚服役。那个国家从法国得到的待遇比突尼斯少得多。法国人建立了对该国及其原住民的统治，但发展了新的通讯、医院和医疗服务。他于1904年6月6日离开阿尔及利亚。

乔尔斯基职业生涯中最重要的一步是他于1905年6月进入陆军地理服务部测地部门，据报道他在那里：-

... 他有敏锐的智力和出色的数学能力，有探究精神和独创性的想法。

德拉姆布雷在1798年春天完成了基线测量，这是他对定义米的贡献的一部分。1882年，法国专家回到了德拉姆布雷的基线，但没有重新测量它，而是倾向于通过三角测量间接检查德拉姆布雷的计算。20世纪第一个十年，巴黎子午线修订后，计划对法国进行新的三角测量。调整网格的问题让地理服务部门的官员非常担忧，他们迫切希望找到一种简单、快速、准确的方法。为了用最小二乘法求解条件方程，Cholesky发明了一个非常巧妙的计算程序，该程序立即被证明非常有用：它现在被称为Cholesky方法，我们将在下面描述它。1905年9月26日，他成为中尉，两年后于1907年5月10日结婚。乔尔斯基和他的妻子有三个孩子；一个儿子和两个女儿。

克里特岛是奥斯曼帝国的一部分已经有200年了，但随着帝国的衰落，欧洲的主要国家打算参与克里特岛的未来。1896年，克里特岛发生了一场反对土耳其人的革命，之后希腊试图取得控制权，列强强行达成了和解。法国人负责西蒂亚区，意大利人负责伊拉佩特拉区，英国人负责伊拉克利安区，俄罗斯人负责雷瑟姆农区，而四人共同负责查尼亚区。法国军队在克里特岛的司令部鲁班斯基上校曾是测地学家，希望在岛上进行测绘工作。1906年3月和4月，测地科长布尔乔亚中校进行了初步侦察。乔尔斯基于1907年11月7日被派往克里特岛，作为执行任务的三名军官之一。他在卡沃西平原测量了基线，并使用天文测量来确定基线南端的精确位置。在他的同事被派去执行其他任务后，乔尔斯基进行了三角测量，并对法国和英国的区域进行了地形测量。他在极其困难的条件下成功地进行了三角测量，因为克里特岛冬季高山上降雪量很大。然而，政治事件阻止了土地调查的进行，乔莱斯基于1908年6月25日离开克里特岛。

1909年3月25日，他被任命为第二步兵团团长，并于1909年9月至1911年9月随部队执行任务。1911年晚些时候，他回到了陆军地理服务局的测地部门，并被分配到阿尔及利亚和突尼斯工作，七年前他已经在那里积累了经验。从1911年10月27日至1912年4月24日，再从1912年10月23日至1913年4月17日，他在阿尔及利亚和突尼斯进行了精密水准测量工作。这项工作旨在为修建铁路线做准备。突尼斯的电网建设于1913-1914年冬季完成，之后有一段时间进行检查和调整。在摩洛哥也必须进行类似的工作，Cholesky设计了一些方法，使工作能够更快地进行，但仍然保持必要的准确性。在测试了他的程序后，这些程序于1912年7月开始在摩洛哥使用，工作于1913年1月完成。

1913年5月25日，乔莱斯基被分配到外交部，并被任命为突尼斯摄政区地形局局长。他开始在那里工作，但在1914年秋天第一次世界大战爆发后，他被调到其他地方工作。1915年1月，他开始组织火炮射击，2月被派往沃斯的一个部队工作。1916年9月，他被派往罗马尼亚执行军事任务，在那里担任罗马尼亚军队地理服务处处长。他于1917年7月晋升，并继续在罗马尼亚服役至1918年2月。

乔莱斯基于1918年8月31日凌晨5点在法国北部的战场上受伤身亡。在他死后，他的一位同事，贝诺伊特司令，在《自然推断法中线性方程组解析法的起源》一书中发表了乔尔斯基的计算方法，以解决一些最小二乘数据拟合问题。应用méthodeála resolution d'un Systemème Definition d'quations lineaires（首席检察官Cholesky）Ⓣ, 发表于1924年的公报géodesique。

Cholesky分解（或Cholesky分解）采用对称正定矩阵AA，并将其写成a=LL'a=LL′

其中，LL是一个下三角矩阵，具有正对角项（有时称为Cholesky三角形），L'L′

是LL的转置。要解Ax=bsolveAx=b，现在需要解LL'x=bLL′

x=b所以把y=L'xy=L′

x表示Ly=bLy=b，解为yy，然后y=L'xy=L′

对x进行求解，得到解。该方法的优点在于，当MM是三角形矩阵时，求解Mx=bMx=b型方程非常简单。

该方法在1924年出版后很少受到关注，但杰克·托德在二战期间将其纳入了伦敦国王学院的分析课程。1948年，福克斯、赫斯基和威尔金森在一篇论文中对该方法进行了分析，而图灵在同一年发表了一篇关于该方法稳定性的论文。

二、乔尔斯基分解（Cholesky decomposition）

矩阵分解 （decomposition, factorization）是将矩阵拆解为数个矩阵的乘积，可分为三角分解、满秩分解、QR分解、Jordan分解和SVD（奇异值）分解等，常见的有三种：
（1）三角分解法（Triangular Factorization）
（2）QR 分解法（QR Factorization）
（3）奇异值分解法（SVD，Singular Value Decomposition）

三角分解法亦称因子分解法，由消元法演变而来的解线性方程组的一类方法。设方程组的矩阵形式为Ax=b，三角分解法就是将系数矩阵A分解为一个下三角矩阵L和一个上三角矩阵U之积：A=LU，然后依次解两个三角形方程组Ly=b和Ux=y，而得到原方程组的解，例如，杜利特尔分解法、乔莱斯基分解法等就是三角分解法。

乔莱斯基分解法（Cholesky decomposition method），亦称平方根法，属于三角分解法之一。

楚列斯基分解（Cholesky decomposition）是1993年公布的数学名词。

乔莱斯基、乔里斯基分解，同。

using System;
using System.Collections;
using System.Collections.Generic;

namespace Legalsoft.Truffer.Algorithm
{
    /// <summary>
    /// 乔莱斯基分解
    /// </summary>
    public static partial class Algorithm_Gallery
    {
        public static int[,] Cholesky_Decomposition(int[,] matrix)
        {
            int n = matrix.GetLength(0);
            int[,] lower = new int[n, n];

            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j <= i; j++)
                {
                    int sum = 0;
                    if (j == i)
                    {
                        for (int k = 0; k < j; k++)
                        {
                            sum += (int)Math.Pow(lower[j, k], 2);
                        }
                        lower[j, j] = (int)Math.Sqrt(matrix[j, j] - sum);
                    }
                    else
                    {
                        for (int k = 0; k < j; k++)
                        {
                            sum += (lower[i, k] * lower[j, k]);
                        }
                        lower[i, j] = (matrix[i, j] - sum) / lower[j, j];
                    }
                }
            }

            return lower;
        }
    }
}