在上一篇文章中，我講述了Logistic回歸的原理以及它的梯度上升法實現(xiàn)。現(xiàn)在來研究Logistic回歸的另一種

實現(xiàn)，即牛頓迭代法。

?

在上篇文章中，我們求出Logistic回歸的似然函數(shù)的偏導數(shù)為

?

??????????????

?

由于是一個多元函數(shù)，變元是，多元函數(shù)求極值問題以前已經(jīng)講過，參考如下文章

?

鏈接：http://blog.csdn.net/acdreamers/article/details/41413787

?

我們知道，極值點的導數(shù)一定均為零，所以一共需要列出個方程，聯(lián)立解出所有的參數(shù)。

當然，這里首先需要用Hessian矩陣來判斷極值的存在性。方程組如下

?

??????????????

?

這一共是個方程，現(xiàn)在的問題變?yōu)槿绾谓膺@個方程組。求Hessian矩陣就得先求二階偏導，即

?

???????????????

?

把Hessian矩陣表示出來，那么有

?

?

所以得到Hessian矩陣，可以看出矩陣是負定的，那么現(xiàn)在我來證明如果是負定的，那

么Hessian矩陣也是負定的。

?

證明：設任意的是n維列向量，因為是負定的，那么為二次型，也是負定的，又因為

?

?????

?

?????所以也是負定的。

?

Hessian矩陣是負定的，也就是說多元函數(shù)存在局部極大值，這符合開始需求的最大似然估計。Hessian

矩陣描述了多元函數(shù)的局部曲率。有了這個Hessian矩陣，我們就可以用牛頓迭代法繼續(xù)進行計算啦！

?

回想一下，對于一元函數(shù)我們是怎樣通過牛頓迭代法求解零點的？ 假設現(xiàn)在要求方程的解，

那么首先選取一個點作為迭代起始點，然后通過下面式子進行迭代，直到達到指定的精度為止。

?

????????????????

?

原理詳見：http://zh.m.wikipedia.org/wiki/%E7%89%9B%E9%A1%BF%E6%B3%95

?

有時候這個起始點的選取很關鍵，因為牛頓迭代法得到的是局部最優(yōu)解，如果函數(shù)只存在一個零點，那么這個

點選取無關重要，但是如果存在多個局部最優(yōu)解，一般是求指定在某個點附近的零點。對于Logistic

回歸問題，Hessian矩陣對于任意數(shù)據(jù)都是負定的，所以說極值點只有一個，初始點選取無關緊要。

?

對于多元函數(shù)求解零點，同樣可以用牛頓迭代法，對于上面的Logistic回歸，可以得到如下迭代式子

?

???????????????

?

其中為Hessian矩陣，而的表示如下

?

??????????????

?

由于Hessian矩陣是對稱負定的，將矩陣提取一個負號出來，得到

?

??????????????

?

然后Hessian矩陣變?yōu)?#xff0c;這樣就是對稱正定的了。那么現(xiàn)在牛頓迭代式變?yōu)?/span>

?

??????????????

?

現(xiàn)在的重點是如何快速并有效計算，即解方程組，通常的做法是直接用高斯消元法求解，

但是這樣做有弊端，弊端有兩個：（1）效率低；（2）數(shù)值穩(wěn)定性差。

?

由于是對稱正定的，可以用Cholesky矩陣分解法來解。Cholesky分解原理如下

?

鏈接：http://blog.csdn.net/acdreamers/article/details/44656847

?

至此，牛頓迭代法求解Logistic回歸的精髓基本講完。現(xiàn)在開始用C++代碼來實現(xiàn)它。

?

代碼：

#include <string.h> #include <fstream> #include <stdio.h> #include <math.h>#include "matrix.h" #define Type double #define Vector vectorusing namespace std;/** 定義數(shù)據(jù)集結構體 */ struct Data {Vector<Type> x;Type y; };/** 預處理數(shù)據(jù)給data */ void PreProcessData(Vector<Data>& data, string path) {string filename = path;ifstream file(filename.c_str());char s[1024];if(file.is_open()){while(file.getline(s, 1024)){Data tmp;Type x1, x2, x3, x4, x5, x6, x7;sscanf(s,"%lf %lf %lf %lf %lf %lf %lf", &x1, &x2, &x3, &x4, &x5, &x6, &x7);tmp.x.push_back(1);tmp.x.push_back(x2);tmp.x.push_back(x3);tmp.x.push_back(x4);tmp.x.push_back(x5);tmp.x.push_back(x6);tmp.y = x7;data.push_back(tmp);}} }void Init(Vector<Data> &data, Vector<Type> &w) {w.clear();data.clear();PreProcessData(data, "TrainData.txt");for(int i = 0; i < data[0].x.size(); i++)w.push_back(0); }Type WX(const Vector<Type>& w, const Data& data) {Type ans = 0;for(int i = 0; i < w.size(); i++)ans += w[i] * data.x[i];return ans; }Type Sigmoid(const Vector<Type>& w, const Data& data) {Type x = WX(w, data);Type ans = exp(x) / (1 + exp(x));return ans; }void PreMatrix(Matrix<Type> &H, Matrix<Type> &U, const Vector<Data> &data, Vector<Type> &w) {int ROWS = data[0].x.size();int COLS = data.size();Matrix<Type> A(COLS, COLS), P(ROWS, COLS), Q(COLS, 1), X(COLS, ROWS);for(int i = 0; i < COLS; i++){Type t = Sigmoid(w, data[i]);A.put(i, i, t *(1 - t));Q.put(i, 0, data[i].y - t);}for(int i = 0; i < ROWS; i++){for(int j = 0; j < COLS; j++)P.put(i, j, data[j].x[i]);}X = P.getTranspose();/** 計算矩陣U和矩陣H的值 */U = P * Q;H = X.getTranspose() * A * X; }Vector<Type> Matrix2Vector(Matrix<Type> &M) {Vector<Type> X;X.clear();int ROWS = M.getRows();for(int i = 0; i < ROWS; i++)X.push_back(M.get(i, 0));return X; }Matrix<Type> Vector2Matrix(Vector<Type> &X) {int ROWS = X.size();Matrix<Type> matrix(ROWS, 1);for(int i = 0; i < ROWS; i++)matrix.put(i, 0, X[i]);return matrix; }/** Cholesky分解得到矩陣L和矩陣D */ void Cholesky(Matrix<Type> &H, Matrix<Type> &L, Matrix<Type> &D) {Type t = 0;int n = H.getRows();for(int k = 0; k < n; k++){for(int i = 0; i < k; i++){t = H.get(i, i) * H.get(k, i) * H.get(k, i);H.put(k, k, H.get(k, k) - t);}for(int j = k + 1; j < n; j++){for(int i = 0; i < k; i++){t = H.get(j, i) * H.get(i, i) * H.get(k, i);H.put(j, k, H.get(j, k) - t);}t = H.get(j, k) / H.get(k, k);H.put(j, k, t);}}for(int i = 0; i < n; i++){D.put(i, i, H.get(i, i));L.put(i, i, 1);for(int j = 0; j < i; j++)L.put(i, j, H.get(i, j));} }/** 回帶求出線性方程組的解 */ void Solve(Matrix<Type> &H, Vector<Type> &X) {int ROWS = H.getRows();int COLS = H.getColumns();Matrix<Type> L(ROWS, COLS), D(ROWS, COLS);Cholesky(H, L, D);int n = ROWS;for(int k = 0; k < n; k++){for(int i = 0; i < k; i++)X[k] -= X[i] * L.get(k, i);X[k] /= L.get(k, k);}L = D * L.getTranspose();for(int k = n - 1; k >= 0; k--){for(int i = k + 1; i < n; i++)X[k] -= X[i] * L.get(k, i);X[k] /= L.get(k, k);} }/** 打印迭代步驟 */ void Display(int cnt, Type error, Vector<Type> w) {cout<<"第"<<cnt<<"次迭代前后的目標差為: "<<error<<endl;cout<<"參數(shù)w為: ";for(int i = 0; i < w.size(); i++)cout<<w[i]<<" ";cout<<endl;cout<<endl; }Type StopFlag(Vector<Type> w1, Vector<Type> w2) {Type ans = 0;int size = w1.size();for(int i = 0; i < size; i++)ans += 0.5 * (w1[i] - w2[i]) * (w1[i] - w2[i]);return ans; }/** 牛頓迭代步驟 */ void NewtonIter(Vector<Data> &data, Vector<Type> &w) {int cnt = 0;Type delta = 0.0001;int ROWS = data[0].x.size();int COLS = data.size();while(1){Matrix<Type> H(ROWS, ROWS), U(ROWS, 1), W(ROWS, 1);PreMatrix(H, U, data, w);Vector<Type> X = Matrix2Vector(U);Solve(H, X);Matrix<Type> x = Vector2Matrix(X);W = Vector2Matrix(w);W += x;Vector<Type> _w = Matrix2Vector(W);Type error = StopFlag(_w, w);w = _w;cnt++;Display(cnt, error, w);if(error < delta) break;} }/** 訓練數(shù)據(jù)得到w數(shù)組，構造分類器 */ void TrainData(Vector<Data> &data, Vector<Type> &w) {Init(data, w);NewtonIter(data, w); }/** 根據(jù)構造好的分類器對數(shù)據(jù)進行分類 */ void Separator(Vector<Type> w) {vector<Data> data;PreProcessData(data, "TestData.txt");cout<<"預測分類結果："<<endl;for(int i = 0; i < data.size(); i++){Type p0 = 0;Type p1 = 0;Type x = WX(w, data[i]);p1 = exp(x) / (1 + exp(x));p0 = 1 - p1;cout<<"實例: ";for(int j = 0; j < data[i].x.size(); j++)cout<<data[i].x[j]<<" ";cout<<"所屬類別為：";if(p1 >= p0) cout<<1<<endl;else cout<<0<<endl;} }int main() {Vector<Type> w;Vector<Data> data;TrainData(data, w);Separator(w);return 0; }

?

上面的牛頓迭代法代碼中用到了矩陣操作，使用頭文件matrix.h，這是一個很優(yōu)秀的矩陣第三方庫，代碼如下

?

代碼：matrix.h

?

/*****************************************************************************/ /* Name: matrix.h */ /* Uses: Class for matrix math functions. */ /* Date: 4/19/2011 */ /* Author: Andrew Que <http://www.DrQue.net/> */ /* Revisions: */ /* 0.1 - 2011/04/19 - QUE - Creation. */ /* 0.5 - 2011/04/24 - QUE - Most functions are complete. */ /* 0.8 - 2011/05/01 - QUE - */ /* = Bug fixes. */ /* + Dot product. */ /* 1.0 - 2011/11/26 - QUE - Release. */ /* */ /* Notes: */ /* This unit implements some very basic matrix functions, which include: */ /* + Addition/subtraction */ /* + Transpose */ /* + Row echelon reduction */ /* + Determinant */ /* + Dot product */ /* + Matrix product */ /* + Scalar product */ /* + Inversion */ /* + LU factorization/decomposition */ /* There isn't much for optimization in this unit as it was designed as */ /* more of a learning experience. */ /* */ /* License: */ /* This program is free software: you can redistribute it and/or modify */ /* it under the terms of the GNU General Public License as published by */ /* the Free Software Foundation, either version 3 of the License, or */ /* (at your option) any later version. */ /* */ /* This program is distributed in the hope that it will be useful, */ /* but WITHOUT ANY WARRANTY; without even the implied warranty of */ /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */ /* GNU General Public License for more details. */ /* */ /* You should have received a copy of the GNU General Public License */ /* along with this program. If not, see <http://www.gnu.org/licenses/>. */ /* */ /* (C) Copyright 2011 by Andrew Que */ /* http://www.DrQue.net/ */ /*****************************************************************************/ #ifndef _MATRIX_H_ #define _MATRIX_H_#include <iostream> #include <cassert> #include <climits> #include <vector>// Class forward for identity matrix. template< class TYPE > class IdentityMatrix;//============================================================================= // Matrix template class // Contains a set of matrix manipulation functions. The template is designed // so that the values of the matrix can be of any type that allows basic // arithmetic. //============================================================================= template< class TYPE = int >class Matrix{protected:// Matrix data.unsigned rows;unsigned columns;// Storage for matrix data.std::vector< std::vector< TYPE > > matrix;// Order sub-index for rows.// Use: matrix[ order[ row ] ][ column ].unsigned * order;//-------------------------------------------------------------// Return the number of leading zeros in the given row.//-------------------------------------------------------------unsigned getLeadingZeros(// Row to countunsigned row) const{TYPE const ZERO = static_cast< TYPE >( 0 );unsigned column = 0;while ( ZERO == matrix[ row ][ column ] )++column;return column;}//-------------------------------------------------------------// Reorder the matrix so the rows with the most zeros are at// the end, and those with the least at the beginning.//// NOTE: The matrix data itself is not manipulated, just the// 'order' sub-indexes.//-------------------------------------------------------------void reorder(){unsigned * zeros = new unsigned[ rows ];for ( unsigned row = 0; row < rows; ++row ){order[ row ] = row;zeros[ row ] = getLeadingZeros( row );}for ( unsigned row = 0; row < (rows-1); ++row ){unsigned swapRow = row;for ( unsigned subRow = row + 1; subRow < rows; ++subRow ){if ( zeros[ order[ subRow ] ] < zeros[ order[ swapRow ] ] )swapRow = subRow;}unsigned hold = order[ row ];order[ row ] = order[ swapRow ];order[ swapRow ] = hold;}delete zeros;}//-------------------------------------------------------------// Divide a row by given value. An elementary row operation.//-------------------------------------------------------------void divideRow(// Row to divide.unsigned row,// Divisor.TYPE const & divisor){for ( unsigned column = 0; column < columns; ++column )matrix[ row ][ column ] /= divisor;}//-------------------------------------------------------------// Modify a row by adding a scaled row. An elementary row// operation.//-------------------------------------------------------------void rowOperation(unsigned row,unsigned addRow,TYPE const & scale){for ( unsigned column = 0; column < columns; ++column )matrix[ row ][ column ] += matrix[ addRow ][ column ] * scale;}//-------------------------------------------------------------// Allocate memory for matrix data.//-------------------------------------------------------------void allocate(unsigned rowNumber,unsigned columnNumber){// Allocate order integers.order = new unsigned[ rowNumber ];// Setup matrix sizes.matrix.resize( rowNumber );for ( unsigned row = 0; row < rowNumber; ++row )matrix[ row ].resize( columnNumber );}//-------------------------------------------------------------// Free memory used for matrix data.//-------------------------------------------------------------void deallocate(unsigned rowNumber,unsigned columnNumber){// Free memory used for storing order (if there is any).if ( 0 != rowNumber )delete[] order;}public:// Used for matrix concatenation.typedef enum{TO_RIGHT,TO_BOTTOM} Position;//-------------------------------------------------------------// Return the number of rows in this matrix.//-------------------------------------------------------------unsigned getRows() const{return rows;}//-------------------------------------------------------------// Return the number of columns in this matrix.//-------------------------------------------------------------unsigned getColumns() const{return columns;}//-------------------------------------------------------------// Get an element of the matrix.//-------------------------------------------------------------TYPE get(unsigned row, // Which row.unsigned column // Which column.) const{assert( row < rows );assert( column < columns );return matrix[ row ][ column ];}//-------------------------------------------------------------// Proform LU decomposition.// This will create matrices L and U such that A=LxU//-------------------------------------------------------------void LU_Decomposition(Matrix & upper,Matrix & lower) const{assert( rows == columns );TYPE const ZERO = static_cast< TYPE >( 0 );upper = *this;lower = *this;for ( unsigned row = 0; row < rows; ++row )for ( unsigned column = 0; column < columns; ++column )lower.matrix[ row ][ column ] = ZERO;for ( unsigned row = 0; row < rows; ++row ){TYPE value = upper.matrix[ row ][ row ];if ( ZERO != value ){upper.divideRow( row, value );lower.matrix[ row ][ row ] = value;}for ( unsigned subRow = row + 1; subRow < rows; ++subRow ){TYPE value = upper.matrix[ subRow ][ row ];upper.rowOperation( subRow, row, -value );lower.matrix[ subRow ][ row ] = value;}}}//-------------------------------------------------------------// Set an element in the matrix.//-------------------------------------------------------------void put(unsigned row,unsigned column,TYPE const & value){assert( row < rows );assert( column < columns );matrix[ row ][ column ] = value;}//-------------------------------------------------------------// Return part of the matrix.// NOTE: The end points are the last elements copied. They can// be equal to the first element when wanting just a single row// or column. However, the span of the total matrix is// ( 0, rows - 1, 0, columns - 1 ).//-------------------------------------------------------------Matrix getSubMatrix(unsigned startRow,unsigned endRow,unsigned startColumn,unsigned endColumn,unsigned const * newOrder = NULL){Matrix subMatrix( endRow - startRow + 1, endColumn - startColumn + 1 );for ( unsigned row = startRow; row <= endRow; ++row ){unsigned subRow;if ( NULL == newOrder )subRow = row;elsesubRow = newOrder[ row ];for ( unsigned column = startColumn; column <= endColumn; ++column )subMatrix.matrix[ row - startRow ][ column - startColumn ] =matrix[ subRow ][ column ];}return subMatrix;}//-------------------------------------------------------------// Return a single column from the matrix.//-------------------------------------------------------------Matrix getColumn(unsigned column){return getSubMatrix( 0, rows - 1, column, column );}//-------------------------------------------------------------// Return a single row from the matrix.//-------------------------------------------------------------Matrix getRow(unsigned row){return getSubMatrix( row, row, 0, columns - 1 );}//-------------------------------------------------------------// Place matrix in reduced row echelon form.//-------------------------------------------------------------void reducedRowEcholon(){TYPE const ZERO = static_cast< TYPE >( 0 );// For each row...for ( unsigned rowIndex = 0; rowIndex < rows; ++rowIndex ){// Reorder the rows.reorder();unsigned row = order[ rowIndex ];// Divide row down so first term is 1.unsigned column = getLeadingZeros( row );TYPE divisor = matrix[ row ][ column ];if ( ZERO != divisor ){divideRow( row, divisor );// Subtract this row from all subsequent rows.for ( unsigned subRowIndex = ( rowIndex + 1 ); subRowIndex < rows; ++subRowIndex ){unsigned subRow = order[ subRowIndex ];if ( ZERO != matrix[ subRow ][ column ] )rowOperation(subRow,row,-matrix[ subRow ][ column ]);}}}// Back substitute all lower rows.for ( unsigned rowIndex = ( rows - 1 ); rowIndex > 0; --rowIndex ){unsigned row = order[ rowIndex ];unsigned column = getLeadingZeros( row );for ( unsigned subRowIndex = 0; subRowIndex < rowIndex; ++subRowIndex ){unsigned subRow = order[ subRowIndex ];rowOperation(subRow,row,-matrix[ subRow ][ column ]);}}} // reducedRowEcholon//-------------------------------------------------------------// Return the determinant of the matrix.// Recursive function.//-------------------------------------------------------------TYPE determinant() const{TYPE result = static_cast< TYPE >( 0 );// Must have a square matrix to even bother.assert( rows == columns );if ( rows > 2 ){int sign = 1;for ( unsigned column = 0; column < columns; ++column ){TYPE subDeterminant;Matrix subMatrix = Matrix( *this, 0, column );subDeterminant = subMatrix.determinant();subDeterminant *= matrix[ 0 ][ column ];if ( sign > 0 )result += subDeterminant;elseresult -= subDeterminant;sign = -sign;}}else{result = ( matrix[ 0 ][ 0 ] * matrix[ 1 ][ 1 ] )- ( matrix[ 0 ][ 1 ] * matrix[ 1 ][ 0 ] );}return result;} // determinant//-------------------------------------------------------------// Calculate a dot product between this and an other matrix.//-------------------------------------------------------------TYPE dotProduct(Matrix const & otherMatrix) const{// Dimentions of each matrix must be the same.assert( rows == otherMatrix.rows );assert( columns == otherMatrix.columns );TYPE result = static_cast< TYPE >( 0 );for ( unsigned row = 0; row < rows; ++row )for ( unsigned column = 0; column < columns; ++column ){result +=matrix[ row ][ column ]* otherMatrix.matrix[ row ][ column ];}return result;} // dotProduct//-------------------------------------------------------------// Return the transpose of the matrix.//-------------------------------------------------------------Matrix const getTranspose() const{Matrix result( columns, rows );// Transpose the matrix by filling the result's rows will// these columns, and vica versa.for ( unsigned row = 0; row < rows; ++row )for ( unsigned column = 0; column < columns; ++column )result.matrix[ column ][ row ] = matrix[ row ][ column ];return result;} // transpose//-------------------------------------------------------------// Transpose the matrix.//-------------------------------------------------------------void transpose(){*this = getTranspose();}//-------------------------------------------------------------// Return inverse matrix.//-------------------------------------------------------------Matrix const getInverse() const{// Concatenate the identity matrix onto this matrix.Matrix inverseMatrix(*this,IdentityMatrix< TYPE >( rows, columns ),TO_RIGHT);// Row reduce this matrix. This will result in the identity// matrix on the left, and the inverse matrix on the right.inverseMatrix.reducedRowEcholon();// Copy the inverse matrix data back to this matrix.Matrix result(inverseMatrix.getSubMatrix(0,rows - 1,columns,columns + columns - 1,inverseMatrix.order));return result;} // invert//-------------------------------------------------------------// Invert this matrix.//-------------------------------------------------------------void invert(){*this = getInverse();} // invert//=======================================================================// Operators.//=======================================================================//-------------------------------------------------------------// Add by an other matrix.//-------------------------------------------------------------Matrix const operator +(Matrix const & otherMatrix) const{assert( otherMatrix.rows == rows );assert( otherMatrix.columns == columns );Matrix result( rows, columns );for ( unsigned row = 0; row < rows; ++row )for ( unsigned column = 0; column < columns; ++column )result.matrix[ row ][ column ] =matrix[ row ][ column ]+ otherMatrix.matrix[ row ][ column ];return result;}//-------------------------------------------------------------// Add self by an other matrix.//-------------------------------------------------------------Matrix const & operator +=(Matrix const & otherMatrix){*this = *this + otherMatrix;return *this;}//-------------------------------------------------------------// Subtract by an other matrix.//-------------------------------------------------------------Matrix const operator -(Matrix const & otherMatrix) const{assert( otherMatrix.rows == rows );assert( otherMatrix.columns == columns );Matrix result( rows, columns );for ( unsigned row = 0; row < rows; ++row )for ( unsigned column = 0; column < columns; ++column )result.matrix[ row ][ column ] =matrix[ row ][ column ]- otherMatrix.matrix[ row ][ column ];return result;}//-------------------------------------------------------------// Subtract self by an other matrix.//-------------------------------------------------------------Matrix const & operator -=(Matrix const & otherMatrix){*this = *this - otherMatrix;return *this;}//-------------------------------------------------------------// Matrix multiplication.//-------------------------------------------------------------Matrix const operator *(Matrix const & otherMatrix) const{TYPE const ZERO = static_cast< TYPE >( 0 );assert( otherMatrix.rows == columns );Matrix result( rows, otherMatrix.columns );for ( unsigned row = 0; row < rows; ++row )for ( unsigned column = 0; column < otherMatrix.columns; ++column ){result.matrix[ row ][ column ] = ZERO;for ( unsigned index = 0; index < columns; ++index )result.matrix[ row ][ column ] +=matrix[ row ][ index ]* otherMatrix.matrix[ index ][ column ];}return result;}//-------------------------------------------------------------// Multiply self by matrix.//-------------------------------------------------------------Matrix const & operator *=(Matrix const & otherMatrix){*this = *this * otherMatrix;return *this;}//-------------------------------------------------------------// Multiply by scalar constant.//-------------------------------------------------------------Matrix const operator *(TYPE const & scalar) const{Matrix result( rows, columns );for ( unsigned row = 0; row < rows; ++row )for ( unsigned column = 0; column < columns; ++column )result.matrix[ row ][ column ] = matrix[ row ][ column ] * scalar;return result;}//-------------------------------------------------------------// Multiply self by scalar constant.//-------------------------------------------------------------Matrix const & operator *=(TYPE const & scalar){*this = *this * scalar;return *this;}//-------------------------------------------------------------// Copy matrix.//-------------------------------------------------------------Matrix & operator =(Matrix const & otherMatrix){if ( this == &otherMatrix )return *this;// Release memory currently in use.deallocate( rows, columns );rows = otherMatrix.rows;columns = otherMatrix.columns;allocate( rows, columns );for ( unsigned row = 0; row < rows; ++row )for ( unsigned column = 0; column < columns; ++column )matrix[ row ][ column ] =otherMatrix.matrix[ row ][ column ];return *this;}//-------------------------------------------------------------// Copy matrix data from array.// Although matrix data is two dimensional, this copy function// assumes the previous row is immediately followed by the next// row's data.//// Example for 3x2 matrix:// int const data[ 3 * 2 ] =// {// 1, 2, 3,// 4, 5, 6// };// Matrix< int > matrix( 3, 2 );// matrix = data;//-------------------------------------------------------------Matrix & operator =(TYPE const * data){unsigned index = 0;for ( unsigned row = 0; row < rows; ++row )for ( unsigned column = 0; column < columns; ++column )matrix[ row ][ column ] = data[ index++ ];return *this;}//-----------------------------------------------------------------------// Return true if this matrix is the same of parameter.//-----------------------------------------------------------------------bool operator ==(Matrix const & value) const{bool isEqual = true;for ( unsigned row = 0; row < rows; ++row )for ( unsigned column = 0; column < columns; ++column )if ( matrix[ row ][ column ] != value.matrix[ row ][ column ] )isEqual = false;return isEqual;}//-----------------------------------------------------------------------// Return true if this matrix is NOT the same of parameter.//-----------------------------------------------------------------------bool operator !=(Matrix const & value) const{return !( *this == value );}//-------------------------------------------------------------// Constructor for empty matrix.// Only useful if matrix is being assigned (i.e. copied) from// somewhere else sometime after construction.//-------------------------------------------------------------Matrix():rows( 0 ),columns( 0 ){allocate( 0, 0 );}//-------------------------------------------------------------// Constructor using rows and columns.//-------------------------------------------------------------Matrix(unsigned rowsParameter,unsigned columnsParameter):rows( rowsParameter ),columns( columnsParameter ){TYPE const ZERO = static_cast< TYPE >( 0 );// Allocate memory for new matrix.allocate( rows, columns );// Fill matrix with zero.for ( unsigned row = 0; row < rows; ++row ){order[ row ] = row;for ( unsigned column = 0; column < columns; ++column )matrix[ row ][ column ] = ZERO;}}//-------------------------------------------------------------// This constructor will allow the creation of a matrix based off// an other matrix. It can copy the matrix entirely, or omitted a// row/column.//-------------------------------------------------------------Matrix(Matrix const & copyMatrix,unsigned omittedRow = INT_MAX,unsigned omittedColumn = INT_MAX){// Start with the number of rows/columns from matrix to be copied.rows = copyMatrix.getRows();columns = copyMatrix.getColumns();// If a row is omitted, then there is one less row.if ( INT_MAX != omittedRow )rows--;// If a column is omitted, then there is one less column.if ( INT_MAX != omittedColumn )columns--;// Allocate memory for new matrix.allocate( rows, columns );unsigned rowIndex = 0;for ( unsigned row = 0; row < rows; ++row ){// If this row is to be skipped...if ( rowIndex == omittedRow )rowIndex++;// Set default order.order[ row ] = row;unsigned columnIndex = 0;for ( unsigned column = 0; column < columns; ++column ){// If this column is to be skipped...if ( columnIndex == omittedColumn )columnIndex++;matrix[ row ][ column ] = copyMatrix.matrix[ rowIndex ][ columnIndex ];columnIndex++;}++rowIndex;}}//-------------------------------------------------------------// Constructor to concatenate two matrices. Concatenation// can be done to the right, or to the bottom.// A = [B | C]//-------------------------------------------------------------Matrix(Matrix const & copyMatrixA,Matrix const & copyMatrixB,Position position = TO_RIGHT){unsigned rowOffset = 0;unsigned columnOffset = 0;if ( TO_RIGHT == position )columnOffset = copyMatrixA.columns;elserowOffset = copyMatrixA.rows;rows = copyMatrixA.rows + rowOffset;columns = copyMatrixA.columns + columnOffset;// Allocate memory for new matrix.allocate( rows, columns );for ( unsigned row = 0; row < copyMatrixA.rows; ++row )for ( unsigned column = 0; column < copyMatrixA.columns; ++column )matrix[ row ][ column ] = copyMatrixA.matrix[ row ][ column ];for ( unsigned row = 0; row < copyMatrixB.rows; ++row )for ( unsigned column = 0; column < copyMatrixB.columns; ++column )matrix[ row + rowOffset ][ column + columnOffset ] =copyMatrixB.matrix[ row ][ column ];}//-------------------------------------------------------------// Destructor.//-------------------------------------------------------------~Matrix(){// Release memory.deallocate( rows, columns );}};//============================================================================= // Class for identity matrix. //============================================================================= template< class TYPE >class IdentityMatrix : public Matrix< TYPE >{public:IdentityMatrix(unsigned rowsParameter,unsigned columnsParameter):Matrix< TYPE >( rowsParameter, columnsParameter ){TYPE const ZERO = static_cast< TYPE >( 0 );TYPE const ONE = static_cast< TYPE >( 1 );for ( unsigned row = 0; row < Matrix< TYPE >::rows; ++row ){for ( unsigned column = 0; column < Matrix< TYPE >::columns; ++column )if ( row == column )Matrix< TYPE >::matrix[ row ][ column ] = ONE;elseMatrix< TYPE >::matrix[ row ][ column ] = ZERO;}}};//----------------------------------------------------------------------------- // Stream operator used to convert matrix class to a string. //----------------------------------------------------------------------------- template< class TYPE >std::ostream & operator<<(// Stream data to place string.std::ostream & stream,// A matrix.Matrix< TYPE > const & matrix){for ( unsigned row = 0; row < matrix.getRows(); ++row ){for ( unsigned column = 0; column < matrix.getColumns(); ++column )stream << "\t" << matrix.get( row , column );stream << std::endl;}return stream;}#endif // _MATRIX_H_

?

總結

以上是生活随笔為你收集整理的Logistic回归与牛顿迭代法的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。