逻辑斯谛回归模型

    逻辑斯谛回归模型是研究因变量为二分类或多分类观察结果与影响因素之间的关系的一种概率型非线性回归模型。逻辑斯谛回归系数通过最大似然估计得到。Logistic函数如下：

         

    式中x为

         

    这里是输入变量的n个特征，然后按照Logistic函数形式求出。

    假设有n个独立变量的向量，设条件概率在x条件下y发生的概率（假设y=1为y发生）。则Logistic函数表示为：

         

    同理，在x条件下y不发生的概率为：

    Logistic回归都是围绕Logistic函数展开的，如何求解是Logistic回归模型的主要问题，采用的最大似然估计来求解这组参数。

    假设有m个观测样本，观测值分别为，设为给定条件下的概率，同理的概率为，得到一个观测值得概率为：

         

    因为各观测样本间相互独立，于是得到似然函数：

         

    对似然函数取对数：

         

    现要求向量使的最大，其中：

         

    要求的最大似然估计，我们需要确定似然函数存在局部极大值。因此，对似然函数求偏导后得：

         

    由多元函数极值的必要条件可知，若多元函数在一点取得极值，且一阶偏导存在，则该点处所有一阶偏导为0。由此，可以得出n+1个方程，如下：

         

    由此方程解出的 不一定是似然函数的极值，需要通过Hessian矩阵来判断得出的解是否为似然函数的极值。

    Hessian矩阵是一个多元函数的二阶偏导构成的方阵，描述了函数的局部曲率。对一个多元函数 ，如果他的二阶偏导都存在，那么Hessian矩阵如下：

         

    通过Hessian矩阵，我们可以判断一点M处极值的三种情况：

1. 如果 是正定矩阵，则临界点M处是一个局部极小值；
2. 如果 是负定矩阵，则临界点M处是一个局部极大值；
3. 如果 是不定矩阵，则临界点M处不是极值。

对于中的n+1个方程，要求Hessian矩阵，先要求似然函数的二阶偏导，即：

         

    则似然函数的Hessian矩阵为

         

    设有矩阵X、A：

         

         

    则似然函数的Hessian矩阵可表示为：

         

    显然，矩阵A是负定的，则可以证明H也是负定的，说明似然函数存在局部极大值。因此，可以使用牛顿迭代法（Newton's Method）来求 。

    对一元函数，使用牛顿迭代法来求零点。假设要求 的解，首先选取一个点 作为迭代起点，通过下面的式子进行迭代，直到达到指定精度为止。

         

    由此，有时起始点选择很关键，如果函数只存在一个零点，那么这个起始点选取就无关重要。对已Logistic回归问题，Hessian矩阵对于任意数据都是负定的，所以说极值点只有一个，初始点的选取无关紧要。

    因此，对于上述Logistic回归的似然函数， 令：

         

    则由可以得到如下的迭代式子：

         

    由于Hessian矩阵是负定的，将矩阵A提取一个负号，得：

         

    然后Hessian矩阵变为

         

    这样，Hessian矩阵就是对称正定的了。那么牛顿迭代式变为：

         

    现在，关键是如何快速并有效的计算，即解方程组 。由于 是对称正定的，可以使用Cholesky矩阵分解法来解。

    若 对称正定，则存在一个对角元为正数的下三角矩阵，使得 成立。对于，可以通过以下步骤求解：

1. 求 的Cholesky分解，得到
2. 求解 ，得到
3. 求解 ，得到

现在的关键问题是对 进行Cholesky分解。假设：

         

    通过比较两边的关系，首先由

         

再由

         

    这样，得到了矩阵 的第一列元素。假设，已经算出了 的前 列元素，通过

         

    可以得出

         

    进一步由

         

    最终：

         

    这样便通过 的前 列求出了第 列，一直递推下去即可求出 。这种方法称为平方根法。

    利用上述方法需要进行开方，这有可能损失精度和增加运算量，为了避免开方，将Cholesky分解进行改进，即：

         

    其中： 是单位下三角矩阵， 为对角均为正数的对角矩阵。把这一分解叫分解。设：

         

    则对于，求解步骤变为：

1. 求 的分解，得到
2. 求解 ，得到
3. 求解 ，得到

对比两边元素，可以得到：

         

    由此可以确定 和 的公式如下：

         

         

         

牛顿迭代法

1. using System;
2. using System.Collections.Generic;
3. using System.Text;
4. using Model.MatrixDecomposition;
5. using Model.Matrix;
6. using System.IO;
7.  
8. namespace Model.NewtonRaphson
9. {
10.     internal class NewtonRaphsonIterate
11.     {
12.         Cholesky_LDL_Decomposition decompositionMatrixByLDL = new Cholesky_LDL_Decomposition();
13.  
14.         private double[,] matrix_Hessian = null;
15.         private double[,] matrix_X = null;
16.         private double[,] matrix_A = null;
17.  
18.  
19.         private double[,] vector_HU = null;
20.         private double[,] vector_U = null;
21.         private double[,] vector_Y = null;
22.         private double[,] vector_Omega = null;
23.         private double[,] vector_PI = null;
24.         private double[,] old_VectorOmega = null;
25.  
26.         #region 属性
27.         public double[,] MatrixL
28.         {
29.             get
30.             {
31.                 return decompositionMatrixByLDL.MatrixL;
32.             }
33.         }
34.  
35.  
36.         public double[,] MatrixD
37.         {
38.             get
39.             {
40.                 return decompositionMatrixByLDL.MatrixD;
41.             }
42.         }
43.  
44.         public double[,] Hessian
45.         {
46.             get
47.             {
48.                 return this.matrix_Hessian;
49.             }
50.         }
51.         #endregion
52.  
53.         /// <summary>
54.         /// 执行牛顿迭代法
55.         /// </summary>
56.         /// <param name="matrix_X">输入矩阵</param>
57.         /// <param name="vector_Y">输出向量</param>
58.         /// <param name="error_Threashold">迭代停止阈值</param>
59.         /// <param name="omega">计算完成的参数</param>
60.         internal void Run(double[,] matrix_X, double[,] vector_Y, double error_Threashold, ref double[,] omega)
61.         {
62.             double error = 1;
63.             this.matrix_X = matrix_X;
64.             this.vector_Y = vector_Y;
65.             this.vector_Omega = omega;
66.             InitializeParameter(matrix_X.GetLength(0));
67.             int i = 0;
68.             while (error > error_Threashold)
69.             {
70.                 i++;
71.                 error = 0;
72.                 old_VectorOmega = (double[,])vector_Omega.Clone();
73.                 GetMatrixAAndVectorPI();
74.                 matrix_Hessian = MatrixOperation.MultipleMatrix(
75.                     MatrixOperation.MatrixMultiDiagMatrix(
76.                     MatrixOperation.TransportMatrix(matrix_X), matrix_A),
77.                     matrix_X);
78.                 GetMatrixU();
79.                 decompositionMatrixByLDL.Cholesky((double[,])matrix_Hessian.Clone(), vector_U, ref vector_HU);
80.                 vector_Omega = MatrixOperation.AddMatrix(vector_Omega, vector_HU);
81.                 GetIterationError(ref error);
82.                 //TODO:迭代计算
83.             }
84.             omega = (double[,])vector_Omega.Clone();
85.         }
86.  
87.         private void InitializeParameter(int rowNumber)
88.         {
89.             matrix_A = new double[rowNumber, 1];
90.             vector_PI = new double[rowNumber, 1];
91.         }
92.  
93.         /// <summary>
94.         /// 获取H=X^TAX的A以及PI(Xi)
95.         /// </summary>
96.         private void GetMatrixAAndVectorPI()
97.         {
98.             for (int i = 0; i < matrix_X.GetLength(0); i++)
99.             {
100.                 double temp_A = 0;
101.                 //matrix_A[i, 0] += vector_Omega[0, 0];
102.                 for (int j = 0; j < matrix_X.GetLength(1); j++)
103.                 {
104.                     temp_A += matrix_X[i, j] * vector_Omega[j, 0];
105.                 }//for2
106.                 matrix_A[i, 0] = (1.0) / (1 + Math.Exp(-temp_A));
107.                 vector_PI[i, 0] = matrix_A[i, 0];
108.                 matrix_A[i, 0] *= (1 - matrix_A[i, 0]);
109.  
110.             }//for1
111.         }
112.  
113.         //计算HU中的U
114.         //U=X^T(Y-PI)
115.         private void GetMatrixU()
116.         {
117.             vector_U = MatrixOperation.MultipleMatrix(MatrixOperation.TransportMatrix(matrix_X),
118.                 MatrixOperation.SubtracteMatrix(vector_Y, vector_PI));
119.         }
120.  
121.         /// <summary>
122.         /// 计算每次迭代完成后的误差
123.         /// </summary>
124.         /// <param name="error">误差</param>
125.         private void GetIterationError(ref double error)
126.         {
127.             for (int i = 0; i < vector_Omega.GetLength(0); i++)
128.             {
129.                 error += Math.Abs(vector_Omega[i, 0] - old_VectorOmega[i, 0]);
130.             }
131.         }
132.     }
133. }

Cholesky分解

1. using System;
2. using System.Collections.Generic;
3. using System.Text;
4. using Model.Matrix;
5.  
6. namespace Model.MatrixDecomposition
7. {
8.     internal class Cholesky_LDL_Decomposition
9.     {
10.         private double[,] matrix_L = null;
11.         private double[,] matrix_D = null;
12.  
13.         public double[,] MatrixL
14.         {
15.             get
16.             {
17.                 return this.matrix_L;
18.             }
19.         }
20.  
21.  
22.         public double[,] MatrixD
23.         {
24.             get
25.             {
26.                 return this.matrix_D;
27.             }
28.         }
29.         private int matrix_Dimension = 0;
30.  
31.         #region Decomposition Matrix
32.         /// <summary>
33.         /// 方程AX=B
34.         /// 利用Cholesky LDL^T分解法，求方程的解
35.         /// </summary>
36.         /// <param name="m_A">系数矩阵A</param>
37.         /// <param name="v_B">列向量，B</param>
38.         /// <param name="v_X">求得方程的解</param>
39.         public void Cholesky(double[,] m_A, double[,] v_B, ref double[,] v_X)
40.         {
41.             this.matrix_Dimension = m_A.GetLength(0);
42.             matrix_L = new double[matrix_Dimension, matrix_Dimension];
43.             matrix_D = new double[matrix_Dimension, matrix_Dimension];
44.             //为了计算方便，将L的严格下三角存储在A的对应位置上，
45.             //而将D的对角元素存储A的对应对角位置上
46.             //double[,] q = (double[,])m_A.Clone();
47.             for (int i = 0; i < matrix_Dimension; i++)
48.             {
49.                 for (int j = 0; j < i; j++)
50.                 {
51.                     m_A[i, i] -= m_A[j, j] * m_A[i, j] * m_A[i, j];
52.                 }//for
53.                 for (int k = i + 1; k < matrix_Dimension; k++)
54.                 {
55.                     for (int n = 0; n < i; n++)
56.                     {
57.                         m_A[k, i] -= m_A[k, n] * m_A[n, n] * m_A[i, n];
58.                     }//for
59.                     m_A[k, i] /= m_A[i, i];
60.                 }//for
61.             }//for
62.  
63.             this.GetLDMatrix(m_A);
64.             this.SolveEquation(v_B,ref v_X);
65.             //double[,] resut = MatrixOperation.MultipleMatrix(MatrixOperation.MultipleMatrix(MatrixL, MatrixD), MatrixOperation.TransportMatrix(MatrixL));
66.         }
67.  
68.         /// <summary>
69.         /// 将L和D矩阵分别赋值
70.         /// </summary>
71.         /// <param name="m_A">经过Cholesky分解后的矩阵</param>
72.         private void GetLDMatrix(double[,] m_A)
73.         {
74.             for (int i = 0; i < matrix_Dimension; i++)
75.             {
76.                 matrix_L[i, i] = 1;
77.                 matrix_D[i, i] = m_A[i, i];
78.                 for (int j = 0; j < i; j++)
79.                 {
80.                     matrix_L[i, j] = m_A[i, j];
81.                 }
82.             }
83.         }
84.         #endregion End Decomposition Matrix
85.  
86.         #region Solve Equation
87.         /// <summary>
88.         /// 求解LY=B => Y
89.         /// DL^TX=Y => X
90.         /// </summary>
91.         /// <param name="v_B">方程AX=B的列向量B</param>
92.         /// <param name="v_X">求解结果</param>
93.         private void SolveEquation(double[,] v_B, ref double[,] v_X)
94.         {
95.             v_X =(double[,])v_B.Clone();
96.             //求解LY=B=>Y
97.             for (int i = 0; i < matrix_Dimension; i++)
98.             {
99.                 for (int j = 0; j < i; j++)
100.                 {
101.                     v_X[i,0] -= v_X[j,0] * matrix_L[i, j];
102.                 }
103.                 v_X[i,0] /= matrix_L[i, i];
104.             }
105.  
106.             //计算DL^T
107.             double[,] dMultiLT = MatrixOperation.MultipleMatrix(matrix_D,
108.                 MatrixOperation.TransportMatrix(matrix_L));
109.             //求解DL^TX=Y => X
110.             for (int i = matrix_Dimension - 1; i >= 0; i--)
111.             {
112.                 for (int j = i + 1; j < matrix_Dimension; j++)
113.                 {
114.                     v_X[i,0] -= v_X[j,0] * dMultiLT[i, j];
115.                 }
116.                 v_X[i,0] /= dMultiLT[i, i];
117.             }
118.         }//end Method
119.         #endregion
120.     }
121. }

MatrixOperation是一个有关矩阵加、减、乘以及特殊矩阵求逆的一个类。