快速用梯度下降法實現一個Logistic Regression 分類器

前陣子聽說一個面試題：你實現一個logistic Regression需要多少分鍾？搞數據挖掘的人都會覺得實現這個簡單的分類器分分鍾就搞定了吧？

因為我做數據挖掘的時候，從來都是順手用用工具的，尤其是微軟內部的TLC相當強大，各種機器學習的算法都有，於是自從離開學校后就沒有自己實現過這些基礎的算法。當有一天心血來潮自己實現一個logistic regression的時候，我會說用了3個小時么？。。。羞羞

---------------------------------------------------前言結束----------------------------------------------

 

當然logistic regression的淵源還是有點深的，想復習理論知識的話可以去http://en.wikipedia.org/wiki/Logistic_regression ， 我這里就直接講實現啦。 

首先要了解一個logistic function

    

這個函數的圖像是這個樣子的：

 

而我們要實現的logistic regression model，就是要去學習出一組權值w：

 x 指feature構成的向量。 這個向量w就可以將每個instance映射到一個實數了。

假如我們要出里的是2分類問題，那么問題就被描述為學習出一組w，使得h(正樣本)趨近於1， h(負樣本)趨近於0.

現在就變成了一個最優化問題，我們要讓誤差最小化。 現在問題來了，怎么定義誤差函數呢？

首先想到的是L2型損失函數啦，於是啪啪啪寫上了

。

很久沒有復習logistic regression的人最容易犯錯的就是在這了。正確的寫法是：

，

然后對它求偏導數得到梯度下降法的迭代更新方程：

。

於是你會發現這個迭代方程是和線性回歸的是一樣的！

理清了過程時候，代碼就變得異常簡單了：

 public class LogisticRegression
    {
        private int _maxIteration = 1000;
        private double _stepSize = 0.000005;
        //private double _stepSize = 0.1;
        private double _lambda = 0.1;
        private double decay = 0.95;

        public int dim;
        public double[] theta;

        public LogisticRegression(int dim)
        {
            this.dim = dim;
        }

        public LogisticRegression(int dim, double stepSize)
            : this(dim)
        {
            this._stepSize = stepSize;
        }

        public void Train(Instance[] instances)
        {
            Initialize();

            int instCnt = instances.Length;
            double[] dev =new double[this.dim];
            for (int t = 0; t < this._maxIteration; t++)
            {
                double cost = 0;
                for (int i = 0; i < instCnt; i++)
                {
                    double h_x = MathLib.Logistic(MathLib.VectorInnerProd(instances[i].featureValues, this.theta));
                    // calculate cost function
                    cost += instances[i].label * Math.Log(h_x) + (1 - instances[i].label) * Math.Log(1 - h_x);
                }
                cost *= -1.0 / instCnt; 
                Console.WriteLine("{0},{1}", t, cost);

               
                for (int i = 0; i < instCnt; i++)
                {
                    ResetArray(dev);
                    double h_x = MathLib.Logistic(MathLib.VectorInnerProd(instances[i].featureValues, this.theta));
                    double error =   h_x- instances[i].label ;
                    for (int j = 0; j < this.dim; j++)
                    {
                        dev[j] += error*instances[i].featureValues[j] + 2*dev[j]*this._lambda;
                        this.theta[j] -= this._stepSize * dev[j] ;
                        //BoundaryLimiting(ref this.theta[j], 0, 1);
                    }
                }
                //this._stepSize *= decay;
                //if (this._stepSize > 0.000001)
                //{
                //    this._stepSize = 0.000001;
                //}
            }
        }
         
        private void BoundaryLimiting(ref double alpha, double lowerbound, double upperbound)
        {
            if (alpha < lowerbound)
            {
                alpha = lowerbound;
            }
            else if (alpha > upperbound)
            {
                alpha = upperbound;
            }
        }
 
        public double[] Predict(Instance[] instances)
        {
            double[] results = new double[instances.Length];
            for (int i = 0; i < results.Length; i++)
            {
                results[i] = MathLib.Logistic(MathLib.VectorInnerProd(instances[i].featureValues, this.theta));
            }
            return results;
        }

        private void ResetArray(double[] dev)
        {
            for (int i = 0; i < dev.Length; i++)
            {
                dev[i] = 0;
            }
        }

        private void Initialize()
        {
            Random ran = new Random(DateTime.Now.Second);

            this.theta = new double[this.dim];
            for (int i = 0; i < this.dim; i++)
            {
                this.theta[i] = ran.NextDouble() * 0 ; // initialize theta with a small value
            }
        }


        public static void Test()
        {
            LogisticRegression lr = new LogisticRegression(3);

            List<Instance> instances = new List<Instance>();
            using (StreamReader rd = new StreamReader(@"D:\\local exp\\data.csv"))
            {
                string content = rd.ReadLine();
                while ((content = rd.ReadLine()) != null)
                {
                    instances.Add(Instance.ParseInstance(content,','));
                }
            }

           // MinMaxNormalize(instances); 

            lr.Train(instances.ToArray()); 
            
        }

        private static void MinMaxNormalize(List<Instance> instances)
        {
            int dim = instances[0].dim;
            double[] min = new double[dim];
            double[] max = new double[dim];

            int instCnt = instances.Count;
            for (int i = 0; i < instCnt; i++)
            {
                for (int j = 0; j < dim; j++)
                {
                    if (i == 0 || instances[i].featureValues[j] < min[j])
                    {
                        min[j] = instances[i].featureValues[j];
                    }
                    if (i == 0 || instances[i].featureValues[j] > max[j])
                    {
                        max[j] = instances[i].featureValues[j];
                    }
                }
            }


            for (int j = 0; j < dim; j++)
            {
                double gap = max[j] - min[j];
                if (gap <= 0)
                {
                    continue;
                }
                for (int i = 0; i < instCnt; i++)
                {
                    instances[i].featureValues[j] = (instances[i].featureValues[j] - min[j]) / gap;
                }
            }
             
        }
    }

前面提到說我花了3個小時，其中很大一部分原因是在debug算法為啥沒有收斂。這里有個很重要的步驟是把feature規范化到[0,1] 。 如果不normalize的話，參數調起來比較麻煩，loss function也經常蹦到NaN去了。

以下是對比normalize和不加normalization的收斂曲線圖：

我用的實現數據可以在 http://pingax.com/wp-content/uploads/2013/12/data.csv  下載到。 它是一個2維的數據， 分布如下：