Logistic回帰とニュートン反復法
39310 ワード
前の記事では、Logistic回帰の原理とその勾配上昇法の実現について述べました.今からLogistic回帰のもう一つを研究します.
実現、すなわちニュートン反復法.
前の文章の中で、Logisticの回帰する尤然関数の偏向導数を求めます.
多元関数なので、元の変化は、多元関数の極値を求める問題です.以前にも述べましたが、次のような文章を参考にしてください.
リンク:http://blog.csdn.net/acdreamers/article/details/41413787
極値点の導関数は必ずゼロであることを知っています.だから、全部で方程式を並べて、連立ですべてのパラメータを解く必要があります.
もちろん、ここではまずHessian行列を用いて極値の存在性を判断する必要がある.方程式グループは以下の通りです
これは全部で方程式です.今の問題はどのようにこの方程式グループを解くかになります.Hessian行列を求めるなら、まず二次偏向ガイドを求めなければなりません.
Hessian行列を表示します.
だから、Hessian行列を得て、行列が負の定であることが分かります.今は負の定であれば、それを証明します.
このHessian行列も負に決められています.
証明:任意のものをn維列ベクトルとします.負のものですから、二次型であり、負のものです.
だから负けです.
Hessian行列は負であり、つまり多元関数には局所的な極大値があり、これは開始ニーズの最大尤度推定に適合する.Hessian
行列は多重関数の局所曲率を記述した.このHessian行列があれば、ニュートン反復法で計算を続けられます.
一要素関数についてはどのようにニュートン反復法で零点を解いていますか?方程式の解が現在要求されていると仮定して、
まず、反復開始点として点を選択し、指定された精度に達するまで、次の式で反復します.
原理詳細は以下の通りですhttp://zh.m.wikipedia.org/wiki/%E7%89%9B%E9%A1%BF%E6%B3%95
この開始点の選択は、ニュートン反復法によって得られた部分的な最適解であり、関数がゼロしか存在しない場合、これは重要である.
点の選択は重要ではないが、複数の局所最適解がある場合は、一般的にある点の近くに指定された零点を求める.Logisticに対して
回帰問題は,Hessian行列は任意のデータに対して負であるので,極値点は一つしかなく,初期点選択は無関係である.
多元関数について零点を解いても、ニュートン反復法を使って、上のLogisticに対して回帰して、次のような反復式が得られます.
この中でHessian行列として、次のように表されています.
Hessian行列は対称負であるため、マトリックスから負の符号を抽出して得られます.
そしてHessian行列が対称的に決まりました.今はニュートンの反復式が
今の重点はどのように迅速かつ効果的に計算するかである.すなわち解方程式群は、一般的な方法はガウスデデメタ解消法で直接的に解いている.
しかし、このようにして弊害があり、弊害が二つあります.(1)効率が低いです.(2)数値安定性が悪い.
対称正定形であるので,Cholesky行列分解法で解くことができる.Cholesky分解の原理は以下の通りです.
リンク:http://blog.csdn.net/acdreamers/article/details/44656847
これでニュートン反復法はLogistic回帰の精髄を解き基本的に説明し終わった.これからC++コードで実現します.
コード:
上のニュートン反復法コードにマトリックス操作が使用されています.先頭ファイルを使用してmarix.hを使います.これは優れたマトリックスサードパーティライブラリです.コードは下記の通りです.
コード:matrix.h
実現、すなわちニュートン反復法.
前の文章の中で、Logisticの回帰する尤然関数の偏向導数を求めます.
多元関数なので、元の変化は、多元関数の極値を求める問題です.以前にも述べましたが、次のような文章を参考にしてください.
リンク:http://blog.csdn.net/acdreamers/article/details/41413787
極値点の導関数は必ずゼロであることを知っています.だから、全部で方程式を並べて、連立ですべてのパラメータを解く必要があります.
もちろん、ここではまずHessian行列を用いて極値の存在性を判断する必要がある.方程式グループは以下の通りです
これは全部で方程式です.今の問題はどのようにこの方程式グループを解くかになります.Hessian行列を求めるなら、まず二次偏向ガイドを求めなければなりません.
Hessian行列を表示します.
だから、Hessian行列を得て、行列が負の定であることが分かります.今は負の定であれば、それを証明します.
このHessian行列も負に決められています.
証明:任意のものをn維列ベクトルとします.負のものですから、二次型であり、負のものです.
だから负けです.
Hessian行列は負であり、つまり多元関数には局所的な極大値があり、これは開始ニーズの最大尤度推定に適合する.Hessian
行列は多重関数の局所曲率を記述した.このHessian行列があれば、ニュートン反復法で計算を続けられます.
一要素関数についてはどのようにニュートン反復法で零点を解いていますか?方程式の解が現在要求されていると仮定して、
まず、反復開始点として点を選択し、指定された精度に達するまで、次の式で反復します.
原理詳細は以下の通りですhttp://zh.m.wikipedia.org/wiki/%E7%89%9B%E9%A1%BF%E6%B3%95
この開始点の選択は、ニュートン反復法によって得られた部分的な最適解であり、関数がゼロしか存在しない場合、これは重要である.
点の選択は重要ではないが、複数の局所最適解がある場合は、一般的にある点の近くに指定された零点を求める.Logisticに対して
回帰問題は,Hessian行列は任意のデータに対して負であるので,極値点は一つしかなく,初期点選択は無関係である.
多元関数について零点を解いても、ニュートン反復法を使って、上のLogisticに対して回帰して、次のような反復式が得られます.
この中でHessian行列として、次のように表されています.
Hessian行列は対称負であるため、マトリックスから負の符号を抽出して得られます.
そしてHessian行列が対称的に決まりました.今はニュートンの反復式が
今の重点はどのように迅速かつ効果的に計算するかである.すなわち解方程式群は、一般的な方法はガウスデデメタ解消法で直接的に解いている.
しかし、このようにして弊害があり、弊害が二つあります.(1)効率が低いです.(2)数値安定性が悪い.
対称正定形であるので,Cholesky行列分解法で解くことができる.Cholesky分解の原理は以下の通りです.
リンク:http://blog.csdn.net/acdreamers/article/details/44656847
これでニュートン反復法はLogistic回帰の精髄を解き基本的に説明し終わった.これからC++コードで実現します.
コード:
#include
#include
#include
#include
#include "matrix.h"
#define Type double
#define Vector vector
using namespace std;
/** */
struct Data
{
Vector x;
Type y;
};
/** data */
void PreProcessData(Vector& 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, Vector &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& 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& w, const Data& data)
{
Type x = WX(w, data);
Type ans = exp(x) / (1 + exp(x));
return ans;
}
void PreMatrix(Matrix &H, Matrix &U, const Vector &data, Vector &w)
{
int ROWS = data[0].x.size();
int COLS = data.size();
Matrix 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 Matrix2Vector(Matrix &M)
{
Vector X;
X.clear();
int ROWS = M.getRows();
for(int i = 0; i < ROWS; i++)
X.push_back(M.get(i, 0));
return X;
}
Matrix Vector2Matrix(Vector &X)
{
int ROWS = X.size();
Matrix matrix(ROWS, 1);
for(int i = 0; i < ROWS; i++)
matrix.put(i, 0, X[i]);
return matrix;
}
/** Cholesky L D */
void Cholesky(Matrix &H, Matrix &L, Matrix &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 &H, Vector &X)
{
int ROWS = H.getRows();
int COLS = H.getColumns();
Matrix 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 w)
{
cout< w1, Vector 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, Vector &w)
{
int cnt = 0;
Type delta = 0.0001;
int ROWS = data[0].x.size();
int COLS = data.size();
while(1)
{
Matrix H(ROWS, ROWS), U(ROWS, 1), W(ROWS, 1);
PreMatrix(H, U, data, w);
Vector X = Matrix2Vector(U);
Solve(H, X);
Matrix x = Vector2Matrix(X);
W = Vector2Matrix(w);
W += x;
Vector _w = Matrix2Vector(W);
Type error = StopFlag(_w, w);
w = _w;
cnt++;
Display(cnt, error, w);
if(error < delta) break;
}
}
/** w ， */
void TrainData(Vector &data, Vector &w)
{
Init(data, w);
NewtonIter(data, w);
}
/** */
void Separator(Vector w)
{
vector data;
PreProcessData(data, "TestData.txt");
cout<= p0) cout<<1< w;
Vector data;
TrainData(data, w);
Separator(w);
return 0;
} 上のニュートン反復法コードにマトリックス操作が使用されています.先頭ファイルを使用してmarix.hを使います.これは優れたマトリックスサードパーティライブラリです.コードは下記の通りです.
コード:matrix.h
/*****************************************************************************/
/* Name: matrix.h */
/* Uses: Class for matrix math functions. */
/* Date: 4/19/2011 */
/* Author: Andrew Que */
/* 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 . */
/* */
/* (C) Copyright 2011 by Andrew Que */
/* http://www.DrQue.net/ */
/*****************************************************************************/
#ifndef _MATRIX_H_
#define _MATRIX_H_
#include
#include
#include
#include
// 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 count
unsigned 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;
else
subRow = 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;
else
result -= 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;
else
rowOffset = 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;
else
Matrix< 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_