FastICA獨立成分 - python實現 - C++實現


FastICA 步驟

1.      對觀測數據 X 進行中心化處理,使樣本的每個屬性均值為0

2.      求出樣本矩陣的協方差矩陣 Cx

3.      用主成分分析得到白化矩陣 W0-1/2UT 對 其中Λ、U分別是Cx的特征值、特征向量

4.      計算正交陣Z=W0*X, Z維數可以根據主成分的占比,小於樣本屬性數

5.      初始化權重矩陣(隨機的)W設迭代次數count、停機誤差Critical

6.      令Wp=E{ Z*g(WTZ) } - E{ g'(WTZ) }*W     其中非線性函數g(x),可取g1(x)=tanh(x)或g2(y)=y*exp(-y2/2)或g3(y)=y3等非線性函數

7.       SZ=Wp*Z 為分離后的信號矩陣


import matplotlib.pyplot as plt  
import numpy as np
from numpy import linalg as LA

C=200 #樣本數
x=np.arange(C)
s1=2*np.sin(0.02*np.pi *x)#正弦信號

a=np.linspace(-2,2,25)
s2= np.array([a,a,a,a,a,a,a,a]).reshape(200,)#鋸齒信號
s3=np.array(20*(5*[2]+5*[-2]))  #方波信號 
s4 =4*(np.random.random([1,C])-0.5).reshape(200,) #隨機信號

# drow origin signal 
ax1 = plt.subplot(411)
ax2 = plt.subplot(412)
ax3 = plt.subplot(413)
ax4 = plt.subplot(414)
ax1.plot(x,s1)
ax2.plot(x,s2)
ax3.plot(x,s3)
ax4.plot(x,s4)
plt.show()

s=np.array([s1,s2,s3,s4])  #合成信號
ran=2*np.random.random([4,4])  #隨機矩陣
mix=ran.dot(s) #混合信號
# drow mix signal 
ax1 = plt.subplot(411)
ax2 = plt.subplot(412)
ax3 = plt.subplot(413)
ax4 = plt.subplot(414)
ax1.plot(x,mix.T[:,0])
ax2.plot(x,mix.T[:,1])
ax3.plot(x,mix.T[:,2])
ax4.plot(x,mix.T[:,3])
plt.show()

#mix=np.array([[1.1,2,3,4,1],
   # [5,6,2,8,7],
  #  [6,2,5,1,3],
   # [7,8,3,3,2]])

Maxcount=10000  #  %最大迭代次數  
Critical=0.00001  #  %判斷是否收斂
R,C=mix.shape

average=np.mean(mix, axis=1) #計算行均值,axis=0,計算每一列的均值

for i in range(R):
    mix[i,:]=mix[i,:]- average[i] #數據標准化,均值為零
Cx=np.cov(mix)  
value,eigvector = np.linalg.eig(Cx)#計算協方差陣的特征值
val=value**(-1/2)*np.eye(R, dtype=float) 
White=np.dot(val ,eigvector.T)  #白化矩陣

Z=np.dot(White,mix) #混合矩陣的主成分Z,Z為正交陣


#W = np.random.random((R,R))# 4x4
W=0.5*np.ones([4,4])#初始化權重矩陣

for n in range(R):
    count=0    
    WP=W[:,n].reshape(R,1) #初始化
    LastWP=np.zeros(R).reshape(R,1) # 列向量;LastWP=zeros(m,1); 
    while LA.norm(WP-LastWP,1)>Critical:
        #print(count," loop :",LA.norm(WP-LastWP,1))
        count=count+1
        LastWP=np.copy(WP)    #  %上次迭代的值  
        gx=np.tanh(LastWP.T.dot(Z))  # 行向量

        for i in range(R):  
            tm1=np.mean( Z[i,:]*gx ) 
            tm2=np.mean(1-gx**2)*LastWP[i] #收斂快
            #tm2=np.mean(gx)*LastWP[i]     #收斂慢
            WP[i]=tm1 - tm2
        #print(" wp :", WP.T ) 
        WPP=np.zeros(R) #一維0向量
        for j in range(n):  
            WPP=WPP+  WP.T.dot(W[:,j])* W[:,j] 
        WP.shape=1,R
        WP=WP-WPP
        WP.shape=R,1
        WP=WP/(LA.norm(WP))
        if(count ==Maxcount):
            print("reach Maxcount,exit loop",LA.norm(WP-LastWP,1))
            break
    print("loop count:",count )
    W[:,n]=WP.reshape(R,)
SZ=W.T.dot(Z)

# plot extract signal
x=np.arange(0,C)
ax1 = plt.subplot(411)
ax2 = plt.subplot(412)
ax3 = plt.subplot(413)
ax4 = plt.subplot(414)
ax1.plot(x, SZ.T[:,0])
ax2.plot(x, SZ.T[:,1])
ax3.plot(x, SZ.T[:,2])
ax4.plot(x, SZ.T[:,3])
plt.show()


結果示意圖:


下面為C++代碼實現,和python有相同的結果。

#include "matrix.h"

/*
數據格式:行數是屬性數,每一列是一個樣本,一下為實例數據
1.1     2       3       4       1
5       6       2       8       7
6       2       5       1       3
7       8       3       3       2
*/

int main()
{
	int R = 4;
	int C = 5;
	int maxcnt= 200;
	Matrix src(R, C);  // 原始數據
	cin >> src;
	src.zeromean(false);

	Matrix Cov = src.T().cov(); //樣本的協方差矩陣
	Matrix eigval = Cov.eig_val();
	cout << "eigval \n" << eigval << endl;
	Matrix eigvect = Cov.eig_vect();
	cout << "eigvect \n" << eigvect << endl;

	eigval = eigval.exponent(-0.5);//每個特征值的 -0.5 次方
	eigval.maxlimit(10000);  ///超過10000,設置為0
	cout << "eigval^-0.5 \n" << eigval;

	Matrix white = eigval* eigvect.T(); //由特征值、特征向量組成白化矩陣
	Matrix Z = white*src;    //正交陣
	cout << "Z:\n" << Z;

	Matrix W(R, R, 0.5);  //初始權值矩陣
	for (size_t i = 0; i < R; i++)
	{
		int cnt = 0;
		Matrix WP = W.getcol(i);  
		Matrix LastWP(R,1,0);
		while ((WP - LastWP).norm1()>0.01)
		{
			cnt++;
			LastWP = WP;  //上次迭代的值
			Matrix gx = (LastWP.T()*Z).mtanh();// 行向量
			for (size_t j = 0; j < R; j++) //更新 WP
			{ 
				double tm1 = (Z.getrow(j).multi(gx) ).mean();
				Matrix one(gx.Row(), gx.Col(), 1.0);
				double tm2 = (one - gx.exponent(2)).mean()*LastWP[j][0];
				WP[j][0] = tm1 - tm2;
			}
			Matrix WPP(R, 1, 0);
			for (size_t k = 0; k < i; k++)
			{
				double tem = (WP.T()*W.getcol(k))[0][0];
				WPP = WPP + tem * W.getcol(k);
			}	
			WP = WP - WPP;
			WP = WP / (WP.norm2());
			if (cnt == maxcnt)
			{
				cout << "reach Maxcount:"<< maxcnt<<",exit loop" << endl;
				break;
			}		
		}
		cout << "loop cnt:" << cnt << endl;
		for (size_t s = 0; s < R; s++)
		{
			W[s][i] = WP[s][0];  
		}
	}
	cout << "W:\n" << W << endl;
	Matrix SZ = W.T()*Z;
	cout << "SZ:\n" << SZ<<endl;
	system("pause");
}

測試數據的結果為:


與python結果對比:



參考文獻:

1,https://en.wikipedia.org/wiki/Independent_component_analysis

2,https://blog.csdn.net/zb1165048017/article/details/48464573


c++引用的矩陣類Matrix.h

#include<iostream>   
#include <fstream>      // std::ifstream
#include <stdlib.h>      
#include <cmath>        
using namespace std;
class Matrix
{
private:
	unsigned  row, col, size;
	double *pmm;//數組指針    
public:
	Matrix(unsigned r, unsigned c) :row(r), col(c)//非方陣構造     
	{
		size = r*c;
		if (size>0)
		{
			pmm = new double[size];
			for (unsigned j = 0; j<size; j++) //init      
			{
				pmm[j] = 0.0;
			}
		}
		else
			pmm = NULL;
	};
	Matrix(unsigned r, unsigned c, double val ) :row(r), col(c)// 賦初值val  
	{
		size = r*c;
		if (size>0)
		{
			pmm = new double[size];
			for (unsigned j = 0; j<size; j++) //init      
			{
				pmm[j] = val;
			}
		}
		else
			pmm = NULL;
	};
	Matrix(unsigned n) :row(n), col(n)//方陣構造      
	{
		size = n*n;
		if (size>0)
		{
			pmm = new double[size];
			for (unsigned j = 0; j<size; j++) //init      
			{
				pmm[j] = 0.0;
			}
		}
		else
			pmm = NULL;
	};
	Matrix(const Matrix &rhs)//拷貝構造   
	{
		row = rhs.row;
		col = rhs.col;
		size = rhs.size;
		pmm = new double[size];
		for (unsigned i = 0; i<size; i++)
			pmm[i] = rhs.pmm[i];
	}

	~Matrix()//析構  
	{
		if (pmm != NULL)
		{
			delete[]pmm;
			pmm = NULL;
		}
	}
 
	Matrix  &operator=(const Matrix&);  //如果類成員有指針必須重寫賦值運算符,必須是成員      
	friend istream &operator>>(istream&, Matrix&);
	
	friend ofstream &operator<<(ofstream &out, Matrix &obj);  // 輸出到文件
	friend ostream &operator<<(ostream&, Matrix&);          // 輸出到屏幕
	friend Matrix  operator+(const Matrix&, const Matrix&);
	friend Matrix  operator-(const Matrix&, const Matrix&);
	friend Matrix  operator*(const Matrix&, const Matrix&);  //矩陣乘法
	friend Matrix  operator*(double, const Matrix&);  //矩陣乘法
	friend Matrix  operator*(const Matrix&, double);  //矩陣乘法

	friend Matrix  operator/(const Matrix&, double);  //矩陣 除以單數

	Matrix multi(const Matrix&); // 對應元素相乘
	Matrix mtanh(); // 對應元素相乘
	unsigned Row()const{ return row; }  
	unsigned Col()const{ return col; }
	Matrix getrow(size_t index); // 返回第index 行,索引從0 算起
	Matrix getcol(size_t index); // 返回第index 列

	Matrix cov(_In_opt_ bool flag = true);   //協方差陣 或者樣本方差   
	double det();   //行列式    
	Matrix solveAb(Matrix &obj);  // b是行向量或者列向量
	Matrix diag();  //返回對角線元素    
	//Matrix asigndiag();  //對角線元素    
	Matrix T()const;   //轉置    
	void sort(bool);//true為從小到大          
	Matrix adjoint();
	Matrix inverse();
	void QR(_Out_ Matrix&, _Out_ Matrix&)const;
	Matrix eig_val(_In_opt_ unsigned _iters = 1000);
	Matrix eig_vect(_In_opt_ unsigned _iters = 1000);

	double norm1();//1范數   
	double norm2();//2范數   
	double mean();// 矩陣均值  
	double*operator[](int i){ return pmm + i*col; }//注意this加括號, (*this)[i][j]   
	void zeromean(_In_opt_  bool flag = true);//默認參數為true計算列
	void normalize(_In_opt_  bool flag = true);//默認參數為true計算列
	Matrix exponent(double x);//每個元素x次冪
	Matrix  eye();//對角陣
	void  maxlimit(double max,double set=0);//對角陣
};
//友元僅僅是指定了訪問權限,不是一般意義的函數聲明,最好在類外再進行一次函數聲明。      
//istream &operator>>(istream &is, Matrix &obj);      
//ostream &operator<<(ostream &is, Matrix &obj);      
//對稱性運算符,如算術,相等,應該是普通非成員函數。      
//Matrix  operator*( const Matrix&,const Matrix& );      
//Matrix  operator+( const Matrix&,const Matrix& );     
//dets遞歸調用    


Matrix Matrix::mtanh() // 對應元素 tanh()
{
	Matrix ret(row, col);
	for (unsigned i = 0; i<ret.size; i++)
	{
		ret.pmm[i] = tanh(pmm[i]);
	}
	return ret;
}
double dets(int n, double *&aa)
{
	if (n == 1)
		return aa[0];
	double *bb = new double[(n - 1)*(n - 1)];//創建n-1階的代數余子式陣bb        
	int mov = 0;//判斷行是否移動       
	double sum = 0.0;//sum為行列式的值      
	for (int arow = 0; arow<n; arow++) // a的行數把矩陣a(nn)賦值到b(n-1)      
	{
		for (int brow = 0; brow<n - 1; brow++)//把aa陣第一列各元素的代數余子式存到bb      
		{
			mov = arow > brow ? 0 : 1; //bb中小於arow的行,同行賦值,等於的錯過,大於的加一      
			for (int j = 0; j<n - 1; j++)  //從aa的第二列賦值到第n列      
			{
				bb[brow*(n - 1) + j] = aa[(brow + mov)*n + j + 1];
			}
		}
		int flag = (arow % 2 == 0 ? 1 : -1);//因為列數為0,所以行數是偶數時候,代數余子式為1.      
		sum += flag* aa[arow*n] * dets(n - 1, bb);//aa第一列各元素與其代數余子式積的和即為行列式    
	}
	delete[]bb;
	return sum;
}

Matrix Matrix::solveAb(Matrix &obj)
{
	Matrix ret(row, 1);
	if (size == 0 || obj.size == 0)
	{
		cout << "solveAb(Matrix &obj):this or obj is null" << endl;
		return ret;
	}
	if (row != obj.size)
	{
		cout << "solveAb(Matrix &obj):the row of two matrix is not equal!" << endl;
		return ret;
	}

	double *Dx = new double[row*row];
	for (int i = 0; i<row; i++)
	{
		for (int j = 0; j<row; j++)
		{
			Dx[i*row + j] = pmm[i*row + j];
		}
	}
	double D = dets(row, Dx);
	if (D == 0)
	{
		cout << "Cramer法則只能計算系數矩陣為滿秩的矩陣" << endl;
		return  ret;
	}


	for (int j = 0; j<row; j++)
	{
		for (int i = 0; i<row; i++)
		{
			for (int j = 0; j<row; j++)
			{
				Dx[i*row + j] = pmm[i*row + j];
			}
		}
		for (int i = 0; i<row; i++)
		{
			Dx[i*row + j] = obj.pmm[i]; //obj賦值給第j列
		}

		//for( int i=0;i<row;i++) //print
		//{
		//	for(int j=0; j<row;j++)
		//	{
		//		cout<< Dx[i*row+j]<<"\t";
		//	} 	
		//	cout<<endl;
		//}
		ret[j][0] = dets(row, Dx) / D;
		 
	}

 
	delete[]Dx;
	return ret;
}

Matrix Matrix::getrow(size_t index)//返回行
{
	Matrix ret(1, col); //一行的返回值
 
	for (unsigned i = 0; i< col; i++)
	{
		 
		 ret[0][i] = pmm[(index) *col + i] ;
	 
	}
	return ret;
}

Matrix Matrix::getcol(size_t index)//返回列
{
	Matrix ret(row, 1); //一列的返回值
 

	for (unsigned i = 0; i< row; i++)
	{

		ret[i][0] = pmm[i *col + index];

	}
	return ret;
}

Matrix Matrix::exponent(double x)//每個元素x次冪
{
	Matrix ret(row, col);
	double a;
	for (unsigned i = 0; i< row; i++)
	{
		for (unsigned j = 0; j < col; j++)
		{
			a=ret[i][j]= pow(pmm[i*col + j],x);
		}
	}
	return ret;
}
void Matrix::maxlimit(double max, double set)//每個元素x次冪
{
 
	for (unsigned i = 0; i< row; i++)
	{
		for (unsigned j = 0; j < col; j++)
		{
			pmm[i*col + j] = pmm[i*col + j]>max ? 0 : pmm[i*col + j];
		}
	}
  
}
Matrix Matrix::eye()//對角陣
{
 
	for (unsigned i = 0; i< row; i++)
	{
		for (unsigned j = 0; j < col; j++)
		{
			if (i == j)
			{
				pmm[i*col + j] = 1.0;
			}
		}
	}
	return *this;
}
void Matrix::zeromean(_In_opt_  bool flag)
{
	if (flag == true) //計算列均值
	{
		double *mean = new double[col];
		for (unsigned j = 0; j < col; j++)
		{
			mean[j] = 0.0;
			for (unsigned i = 0; i < row; i++)
			{
				mean[j] += pmm[i*col + j];
			}
			mean[j] /= row;
		}
		for (unsigned j = 0; j < col; j++)
		{

			for (unsigned i = 0; i < row; i++)
			{
				pmm[i*col + j] -= mean[j];
			}
		}
		delete[]mean;
	}
	else //計算行均值
	{
		double *mean = new double[row];
		for (unsigned i = 0; i< row; i++)
		{
			mean[i] = 0.0;
			for (unsigned j = 0; j < col; j++)
			{
				mean[i] += pmm[i*col + j];
			}
			mean[i] /= col;
		}
		for (unsigned i = 0; i < row; i++)
		{
			for (unsigned j = 0; j < col; j++)
			{
				pmm[i*col + j] -= mean[i];
			}
		}
		delete[]mean;
	}
}

void Matrix::normalize(_In_opt_  bool flag)
{
	if (flag == true) //計算列均值
	{
		double *mean = new double[col];

		for (unsigned j = 0; j < col; j++)
		{
			mean[j] = 0.0;
			for (unsigned i = 0; i < row; i++)
			{
				mean[j] += pmm[i*col + j];
			}
			mean[j] /= row;
		}
		for (unsigned j = 0; j < col; j++)
		{

			for (unsigned i = 0; i < row; i++)
			{
				pmm[i*col + j] -= mean[j];
			}
		}
		///計算標准差
		for (unsigned j = 0; j < col; j++)
		{
			mean[j] = 0;
			for (unsigned i = 0; i < row; i++)
			{
				mean[j] += pow(pmm[i*col + j],2);//列平方和
			}
				mean[j] = sqrt(mean[j] / row); // 開方
		}
		for (unsigned j = 0; j < col; j++)
		{
			for (unsigned i = 0; i < row; i++)
			{
				pmm[i*col + j] /= mean[j];//列平方和
			}
		}
		delete[]mean;
	}
	else //計算行均值
	{
		double *mean = new double[row];
		for (unsigned i = 0; i< row; i++)
		{
			mean[i] = 0.0;
			for (unsigned j = 0; j < col; j++)
			{
				mean[i] += pmm[i*col + j];
			}
			mean[i] /= col;
		}
		for (unsigned i = 0; i < row; i++)
		{
			for (unsigned j = 0; j < col; j++)
			{
				pmm[i*col + j] -= mean[i];
			}
		}
		///計算標准差
		for (unsigned i = 0; i< row; i++)
		{
			mean[i] = 0.0;
			for (unsigned j = 0; j < col; j++)
			{
				mean[i] += pow(pmm[i*col + j], 2);//列平方和
			}
			mean[i] = sqrt(mean[i] / col); // 開方
		}
		for (unsigned i = 0; i < row; i++)
		{
			for (unsigned j = 0; j < col; j++)
			{
				pmm[i*col + j] /= mean[i];
			}
		}
		delete[]mean;
	}
}

double Matrix::det()
{
	if (col == row)
		return dets(row, pmm);
	else
	{
		cout << ("行列不相等無法計算") << endl;
		return 0;
	}
}
/////////////////////////////////////////////////////////////////////      
istream &operator>>(istream &is, Matrix &obj)
{
	for (unsigned i = 0; i<obj.size; i++)
	{
		is >> obj.pmm[i];
	}
	return is;
}

ostream &operator<<(ostream &out, Matrix &obj)
{
	for (unsigned i = 0; i < obj.row; i++) //打印逆矩陣        
	{
		for (unsigned j = 0; j < obj.col; j++)
		{
			out << (obj[i][j]) << "\t";
		}
		out << endl;
	}
	return out;
}
ofstream &operator<<(ofstream &out, Matrix &obj)//打印逆矩陣到文件  
{
	for (unsigned i = 0; i < obj.row; i++)      
	{
		for (unsigned j = 0; j < obj.col; j++)
		{
			out << (obj[i][j]) << "\t";
		}
		out << endl;
	}
	return out;
}
Matrix operator+(const Matrix& lm, const Matrix& rm)
{
	if (lm.col != rm.col || lm.row != rm.row)
	{
		Matrix temp(0, 0);
		temp.pmm = NULL;
		cout << "operator+(): 矩陣shape 不合適,col:"
			<< lm.col << "," << rm.col << ".  row:" << lm.row << ", " << rm.row << endl;
		return temp; //數據不合法時候,返回空矩陣      
	}
	Matrix ret(lm.row, lm.col);
	for (unsigned i = 0; i<ret.size; i++)
	{
		ret.pmm[i] = lm.pmm[i] + rm.pmm[i];
	}
	return ret;
}
Matrix operator-(const Matrix& lm, const Matrix& rm)
{
	if (lm.col != rm.col || lm.row != rm.row)
	{
		Matrix temp(0, 0);
		temp.pmm = NULL;
		cout << "operator-(): 矩陣shape 不合適,col:" 
			<<lm.col<<","<<rm.col<<".  row:"<< lm.row <<", "<< rm.row << endl;

		return temp; //數據不合法時候,返回空矩陣      
	}
	Matrix ret(lm.row, lm.col);
	for (unsigned i = 0; i<ret.size; i++)
	{
		ret.pmm[i] = lm.pmm[i] - rm.pmm[i];
	}
	return ret;
}
Matrix operator*(const Matrix& lm, const Matrix& rm)  //矩陣乘法
{
	if (lm.size == 0 || rm.size == 0 || lm.col != rm.row)
	{
		Matrix temp(0, 0);
		temp.pmm = NULL;
		cout << "operator*(): 矩陣shape 不合適,col:"
			<< lm.col << "," << rm.col << ".  row:" << lm.row << ", " << rm.row << endl;
		return temp; //數據不合法時候,返回空矩陣      
	}
	Matrix ret(lm.row, rm.col);
	for (unsigned i = 0; i<lm.row; i++)
	{
		for (unsigned j = 0; j< rm.col; j++)
		{
			for (unsigned k = 0; k< lm.col; k++)//lm.col == rm.row      
			{
				ret.pmm[i*rm.col + j] += lm.pmm[i*lm.col + k] * rm.pmm[k*rm.col + j];
			}
		}
	}
	return ret;
}
Matrix operator*(double val, const Matrix& rm)  //矩陣乘 單數
{
	Matrix ret(rm.row, rm.col);
	for (unsigned i = 0; i<ret.size; i++)
	{
		ret.pmm[i] = val * rm.pmm[i];
	}
	return ret;
}
Matrix operator*(const Matrix&lm, double val)  //矩陣乘 單數
{
	Matrix ret(lm.row, lm.col);
	for (unsigned i = 0; i<ret.size; i++)
	{
		ret.pmm[i] = val * lm.pmm[i];
	}
	return ret;
}

Matrix operator/(const Matrix&lm, double val)  //矩陣除以 單數
{
	Matrix ret(lm.row, lm.col);
	for (unsigned i = 0; i<ret.size; i++)
	{
		ret.pmm[i] =  lm.pmm[i]/val;
	}
	return ret;
}
Matrix Matrix::multi(const Matrix&rm)// 對應元素相乘
{
	if (col != rm.col || row != rm.row)
	{
		Matrix temp(0, 0);
		temp.pmm = NULL;
		cout << "multi(const Matrix&rm): 矩陣shape 不合適,col:"
			<< col << "," << rm.col << ".  row:" << row << ", " << rm.row << endl;
		return temp; //數據不合法時候,返回空矩陣      
	}
	Matrix ret(row,col);
	for (unsigned i = 0; i<ret.size; i++)
	{
		ret.pmm[i] = pmm[i] * rm.pmm[i];
	}
	return ret;

}

Matrix&  Matrix::operator=(const Matrix& rhs)
{
	if (this != &rhs)
	{
		row = rhs.row;
		col = rhs.col;
		size = rhs.size;
		if (pmm != NULL)
			delete[] pmm;
		pmm = new double[size];
		for (unsigned i = 0; i<size; i++)
		{
			pmm[i] = rhs.pmm[i];
		}
	}
	return *this;
}
//||matrix||_2  求A矩陣的2范數        
double Matrix::norm2()
{
	double norm = 0;
	for (unsigned i = 0; i < size; ++i)
	{
		norm += pmm[i] * pmm[i];
	}
	return (double)sqrt(norm);
}
double Matrix::norm1()
{
	double sum = 0;
	for (unsigned i = 0; i < size; ++i)
	{
		sum += abs(pmm[i]);
	}
	return sum;
}
double Matrix::mean()
{
	double sum = 0;
	for (unsigned i = 0; i < size; ++i)
	{
		sum += (pmm[i]);
	}
	return sum/size;
}




void Matrix::sort(bool flag)
{
	double tem;
	for (unsigned i = 0; i<size; i++)
	{
		for (unsigned j = i + 1; j<size; j++)
		{
			if (flag == true)
			{
				if (pmm[i]>pmm[j])
				{
					tem = pmm[i];
					pmm[i] = pmm[j];
					pmm[j] = tem;
				}
			}
			else
			{
				if (pmm[i]<pmm[j])
				{
					tem = pmm[i];
					pmm[i] = pmm[j];
					pmm[j] = tem;
				}
			}

		}
	}
}
Matrix Matrix::diag()
{
	if (row != col)
	{
		Matrix m(0);
		cout << "diag():row != col" << endl;
		return m;
	}
	Matrix m(row);
	for (unsigned i = 0; i<row; i++)
	{
		m.pmm[i*row + i] = pmm[i*row + i];
	}
	return m;
}
Matrix Matrix::T()const
{
	Matrix tem(col, row);
	for (unsigned i = 0; i<row; i++)
	{
		for (unsigned j = 0; j<col; j++)
		{
			tem[j][i] = pmm[i*col + j];// (*this)[i][j]  
		}
	}
	return tem;
}
void  Matrix::QR(Matrix &Q, Matrix &R) const
{
	//如果A不是一個二維方陣,則提示錯誤,函數計算結束        
	if (row != col)
	{
		printf("ERROE: QR() parameter A is not a square matrix!\n");
		return;
	}
	const unsigned N = row;
	double *a = new double[N];
	double *b = new double[N];

	for (unsigned j = 0; j < N; ++j)  //(Gram-Schmidt) 正交化方法      
	{
		for (unsigned i = 0; i < N; ++i)  //第j列的數據存到a,b      
			a[i] = b[i] = pmm[i * N + j];

		for (unsigned i = 0; i<j; ++i)  //第j列之前的列      
		{
			R.pmm[i * N + j] = 0;  //      
			for (unsigned m = 0; m < N; ++m)
			{
				R.pmm[i * N + j] += a[m] * Q.pmm[m *N + i]; //R[i,j]值為Q第i列與A的j列的內積      
			}
			for (unsigned m = 0; m < N; ++m)
			{
				b[m] -= R.pmm[i * N + j] * Q.pmm[m * N + i]; //      
			}
		}

		double norm = 0;
		for (unsigned i = 0; i < N; ++i)
		{
			norm += b[i] * b[i];
		}
		norm = (double)sqrt(norm);

		R.pmm[j*N + j] = norm; //向量b[]的2范數存到R[j,j]      

		for (unsigned i = 0; i < N; ++i)
		{
			Q.pmm[i * N + j] = b[i] / norm; //Q 陣的第j列為單位化的b[]      
		}
	}
	delete[]a;
	delete[]b;
}
Matrix Matrix::eig_val(_In_opt_ unsigned _iters)
{
	if (size == 0 || row != col)
	{
		cout << "矩陣為空或者非方陣!" << endl;
		Matrix rets(0);
		return rets;
	}
	//if (det() == 0)  
	//{  
	//  cout << "非滿秩矩陣沒法用QR分解計算特征值!" << endl;  
	//  Matrix rets(0);  
	//  return rets;  
	//}  
	const unsigned N = row;
	Matrix matcopy(*this);//備份矩陣      
	Matrix Q(N), R(N);
	/*當迭代次數足夠多時,A 趨於上三角矩陣,上三角矩陣的對角元就是A的全部特征值。*/
	for (unsigned k = 0; k < _iters; ++k)
	{
		//cout<<"this:\n"<<*this<<endl;      
		QR(Q, R);
		*this = R*Q;
		/*  cout<<"Q:\n"<<Q<<endl;
		cout<<"R:\n"<<R<<endl;  */
	}
	Matrix val = diag();
	*this = matcopy;//恢復原始矩陣;    
	return val;
}
Matrix Matrix::eig_vect(_In_opt_ unsigned _iters)
{
	if (size == 0 || row != col)
	{
		cout << "矩陣為空或者非方陣!" << endl;
		Matrix rets(0);
		return rets;
	}
	if (det() == 0)  
	{  
	  cout << "非滿秩矩陣沒法用QR分解計算特征向量!" << endl;  
	  Matrix rets(0);  
	  return rets;  
	}  
	Matrix matcopy(*this);//備份矩陣     
	Matrix eigenValue = eig_val(_iters);
	Matrix ret(row);
	const unsigned NUM = col;
	double eValue;
	double sum, midSum, diag;
	Matrix copym(*this);
	for (unsigned count = 0; count < NUM; ++count)
	{
		//計算特征值為eValue,求解特征向量時的系數矩陣       
		*this = copym;
		eValue = eigenValue[count][count];

		for (unsigned i = 0; i < col; ++i)//A-lambda*I       
		{
			pmm[i * col + i] -= eValue;
		}
		//cout<<*this<<endl;      
		//將 this為階梯型的上三角矩陣        
		for (unsigned i = 0; i < row - 1; ++i)
		{
			diag = pmm[i*col + i];  //提取對角元素    
			for (unsigned j = i; j < col; ++j)
			{
				pmm[i*col + j] /= diag; //【i,i】元素變為1    
			}
			for (unsigned j = i + 1; j<row; ++j)
			{
				diag = pmm[j *  col + i];
				for (unsigned q = i; q < col; ++q)//消去第i+1行的第i個元素    
				{
					pmm[j*col + q] -= diag*pmm[i*col + q];
				}
			}
		}
		//cout<<*this<<endl;      
		//特征向量最后一行元素置為1    
		midSum = ret.pmm[(ret.row - 1) * ret.col + count] = 1;
		for (int m = row - 2; m >= 0; --m)
		{
			sum = 0;
			for (unsigned j = m + 1; j < col; ++j)
			{
				sum += pmm[m *  col + j] * ret.pmm[j * ret.col + count];
			}
			sum = -sum / pmm[m *  col + m];
			midSum += sum * sum;
			ret.pmm[m * ret.col + count] = sum;
		}
		midSum = sqrt(midSum);
		for (unsigned i = 0; i < ret.row; ++i)
		{
			ret.pmm[i * ret.col + count] /= midSum; //每次求出一個列向量      
		}
	}
	*this = matcopy;//恢復原始矩陣;    
	return ret;
}
Matrix Matrix::cov(bool flag)
{
	//row 樣本數,column 變量數    
	if (col == 0)
	{
		Matrix m(0);
		return m;
	}
	double *mean = new double[col]; //均值向量    

	for (unsigned j = 0; j<col; j++) //init    
	{
		mean[j] = 0.0;
	}
	Matrix ret(col);
	for (unsigned j = 0; j<col; j++) //mean    
	{
		for (unsigned i = 0; i<row; i++)
		{
			mean[j] += pmm[i*col + j];
		}
		mean[j] /= row;
	}
	unsigned i, k, j;
	for (i = 0; i<col; i++) //第一個變量    
	{
		for (j = i; j<col; j++) //第二個變量    
		{
			for (k = 0; k<row; k++) //計算    
			{
				ret[i][j] += (pmm[k*col + i] - mean[i])*(pmm[k*col + j] - mean[j]);

			}
			if (flag == true)
			{
				ret[i][j] /= (row-1);
			}
			else
			{
				ret[i][j] /= (row);
			}
			/*temp.Format("cov=%f,column=%d mean(%d)=%f,mean(%d)=%f",cov[i*column+j],row,i,mean[i],j,mean[j]);
			MessageBox(temp);*/
		}
	}
	for (i = 0; i<col; i++) //補全對應面    
	{
		for (j = 0; j<i; j++)
		{
			ret[i][j] = ret[j][i];
		}
	}
	return ret;
}
Matrix Matrix::adjoint()
{
	//調動之前,檢查時候方陣,這里默認為aa為方陣    
	//確定本函數是否修改傳入的數據 :no    
	//調用函數內刪除內存delete []padjoint;    
	if (row != col)
	{
		Matrix ret(0);
		return ret;
	}

	const int n = row;
	Matrix ret(row);
	double *bb = new double[(n - 1)*(n - 1)];//創建n-1階的代數余子式陣bb     

	int pi, pj, q;
	for (int ai = 0; ai<n; ai++) // a的行數把矩陣a(nn)賦值到b(n-1)    
	{
		for (int aj = 0; aj<n; aj++)
		{
			for (int bi = 0; bi<n - 1; bi++)//把元素aa[ai][0]余子式存到bb[][]    
			{
				for (int bj = 0; bj<n - 1; bj++)//把元素aa[ai][0]代數余子式存到bb[][]    
				{
					if (ai>bi)    //ai行的代數余子式是:小於ai的行,aa與bb陣,同行賦值    
						pi = 0;
					else
						pi = 1;     //大於等於ai的行,取aa陣的ai+1行賦值給陣bb的bi行    
					if (aj>bj)    //ai行的代數余子式是:小於ai的行,aa與bb陣,同行賦值    
						pj = 0;
					else
						pj = 1;     //大於等於ai的行,取aa陣的ai+1行賦值給陣bb的bi行    
					bb[bi*(n - 1) + bj] = pmm[(bi + pi)*n + bj + pj];
				}
			}
			if ((ai + aj) % 2 == 0)  q = 1;//因為列數為0,所以行數是偶數時候,代數余子式為-1.    
			else  q = (-1);
			ret.pmm[ai*n + aj] = q*dets(n - 1, bb);  //加符號變為代數余子式    
		}
	}
	delete[]bb;
	return ret;
}
Matrix Matrix::inverse()
{
	double det_aa = det();
	if (det_aa == 0)
	{
		cout << "行列式為0 ,不能計算逆矩陣。" << endl;
		Matrix rets(0);
		return rets;
	}
	Matrix adj = adjoint();
	Matrix ret(row);

	for (unsigned i = 0; i<row; i++) //print    
	{
		for (unsigned j = 0; j<col; j++)
		{
			ret.pmm[i*col + j] = adj.pmm[i*col + j] / det_aa;
		}
	}
	return ret;
}



注意!

本站转载的文章为个人学习借鉴使用,本站对版权不负任何法律责任。如果侵犯了您的隐私权益,请联系我们删除。



 
粤ICP备14056181号  © 2014-2021 ITdaan.com