matrix.h

矩阵运算的模板类

// Conjugate Transpose: return conjugate transposed matrix of given matrix.
template <typename T> Matrix<T> Matrix<T>::ct() const {
	int i,j;
	Matrix<T> B(n, m);
	for (i=1; i<=m; i++) for (j=1; j<=n; j++){
			B.mat[j][i]=conj(mat[i][j]);
	}
	return B;
}


template <typename T> inline Matrix<T> Matrix<T>::i() const {
	return Inv(*this);
}


template <typename T> inline int Matrix<T>::p(string& name) const {
	return (*this).Print(cout,name);
}


template <typename T> inline Matrix<T> Transpose(const Matrix<T>& A) {
	return A.ct();
}


template <typename T> inline Matrix<T> ArrayTranspose(const Matrix<T>& A) {
	return A.ct();
}


// Conjugate: return conjugate transposed matrix of given matrix.
template <typename T> Matrix<T> Conj(const Matrix<T>& A) {
	int i,j;
	Matrix<T> B(A.m, A.n);
	for (i=1; i<=m; i++) for (j=1; j<=n; j++){
			B.mat[i][j]=conj(mat[i][j]);
	}
	return B;
}


// return the real matrix of a complex matrix.
template <typename T> Matrix<double> Matrix<T>::Real() const {
	int i,j;
	Matrix<double> B(m, n);
	for (i=1; i<=m; i++) for (j=1; j<=n; j++){
			B[i][j]=real((*this)[i][j]);
	}
	return B;
}


// return the imaginary matrix of a complex matrix.
template <typename T> Matrix<double> Matrix<T>::Imag() const {
	int i,j;
	Matrix<double> B(m, n);
	for (i=1; i<=m; i++) for (j=1; j<=n; j++){
			B[i][j]=imag((*this)[i][j]);
	}
	return B;
}


template <typename T> Vector<T> Diag(const Matrix<T>& A){
	Vector<T> C(A.m,COL_VECTOR);
	for(int i=1;i<=A.m;i++) C[i]=A.mat[i][i];
	return C;
}


template <typename T> inline Vector<T> Sum(const Matrix<T>& A,int idx=1){
	return A.Sum_(idx);
}


template <typename T> Vector<T> Mean(const Matrix<T>& A,int idx=1){
	assert(idx==1 || idx==2);
	int sampleN=(idx==1)?NumRows(A):NumCols(A);
	return Sum(A,idx)/T(sampleN);
}


template <typename T> Matrix<T> Cov(const Matrix<T>& A){
	Vector<T> mean=Mean(A);
	return (Transpose(A)*A)/(double)(A.m)-OuterProd(mean,mean);
}


// Trace
template <typename T> T Trace(const Matrix<T>& A){
	assert(A.m==A.n);
	T sum=T();
	for(int i=1;i<=A.m;i++) sum+=A[i][i];
	return T(sum);
}


// Array Multiply: A(i,j)*B(i,j)
template <typename T> Matrix<T>  ArrayMultiply(const Matrix<T>& A, const Matrix<T>& B){
	return mtl_times(A,B);
}


// Array Division: A(i,j)/B(i,j)
template <typename T> Matrix<T>  ArrayDivide(const Matrix<T>& A, const Matrix<T>& B){
	return mtl_rdivide(A,B);
}


template <typename T> Matrix<T> Abs(const Matrix<T>& A){
	Matrix<T> C(A.m,A.n);
	for(int i=1;i<=A.m;i++)
		for(int j=1;j<=A.n;j++) C[i][j]=T(Abs(A[i][j]));
	return C;
}


template <typename T> Matrix<T> Sign(const Matrix<T>& A, double x){
	Matrix<T> C(A.m,A.n);
	for(int i=1;i<=A.m;i++)
		for(int j=1;j<=A.n;j++) C[i][j]=T( ( (A[i][j]>0) ? 1 : ((A[i][j]<0) ? (-1) : 0) ) );
	return C;
}


template <typename T> Matrix<T> Pow(const Matrix<T>& A, double x){
	return mtl_power(A,x);
}


template <typename T> Matrix<T> Pow(const Matrix<T>& A, int x){
	return mtl_power(A,x);
}


template <typename T> Matrix<T> Powm(const Matrix<T>& A, double x){
	return mtl_mpower(A,x);
}


template <typename T> Matrix<T> Powm(const Matrix<T>& A, int x){
	return mtl_mpower(A,x);
}


template <typename T> Matrix<T> Exp(const Matrix<T>& A){
	Matrix<T> C(A.m,A.n);
	for(int i=1;i<=A.m;i++)
		for(int j=1;j<=A.n;j++) C[i][j]=exp(A[i][j]);
	return C;
}


template <typename T> Matrix<T> Log(const Matrix<T>& A){
	Matrix<T> C(A.m,A.n);
	for(int i=1;i<=A.m;i++)
		for(int j=1;j<=A.n;j++) C[i][j]=log(A[i][j]);
	return C;
}


template <typename T> Matrix<T> Sqrt(const Matrix<T>& A){
	Matrix<T> C(A.m,A.n);
	for(int i=1;i<=A.m;i++)
		for(int j=1;j<=A.n;j++) C[i][j]=sqrt(A[i][j]);
	return C;
}


template <typename T> void Matrix<T>::FromInt(const Matrix<int>& A){
	Resize(NumRows(A),NumCols(A));
	for(int i=1;i<=m;i++)
		for(int j=1;j<=n;j++) (*this)[i][j]=T(A[i][j]);
}


template <typename T> void Matrix<T>::FromFloat(const Matrix<float>& A){
	Resize(NumRows(A),NumCols(A));
	for(int i=1;i<=m;i++)
		for(int j=1;j<=n;j++) (*this)[i][j]=T(A[i][j]);
}


template <typename T> void Matrix<T>::FromDouble(const Matrix<double>& A){
	Resize(NumRows(A),NumCols(A));
	for(int i=1;i<=m;i++)
		for(int j=1;j<=n;j++) (*this)[i][j]=T(A[i][j]);
}


#ifndef DISABLE_COMPLEX
inline void Matrix<complex<float> >::FromCDouble(const Matrix<complex<double> >& A){
	Resize(NumRows(A),NumCols(A));
	for(int i=1;i<=m;i++)
		for(int j=1;j<=n;j++) (*this)[i][j]=complex<float>(real(A[i][j]),imag(A[i][j]));
}
#endif


template <typename T> Matrix<int> Matrix<T>::ToInt() const {
	Matrix<int> C(m,n);
	for(int i=1;i<=m;i++)
		for(int j=1;j<=n;j++) C[i][j]=(int)(*this)[i][j];
	return C;
}


template <typename T> Matrix<float> Matrix<T>::ToFloat() const {
	Matrix<float> C(m,n);
	for(int i=1;i<=m;i++)
		for(int j=1;j<=n;j++) C[i][j]=(float)(*this)[i][j];
	return C;
}


template <typename T> Matrix<double> Matrix<T>::ToDouble() const {
	Matrix<double> C(m,n);
	for(int i=1;i<=m;i++)
		for(int j=1;j<=n;j++) C[i][j]=(double)(*this)[i][j];
	return C;
}


#ifndef DISABLE_COMPLEX
template <typename T> Matrix<complex<double> > Matrix<T>::ToCDouble() const {
	Matrix<complex<double> > C(m,n);
	for(int i=1;i<=m;i++)
		for(int j=1;j<=n;j++) C[i][j]=complex<double>(real((*this)[i][j]), imag((*this)[i][j]));
	return C;
}
#endif


template <typename T> Matrix<int> Int(const Matrix<T>& A) {
	Matrix<int> C(NumRows(A),NumCols(A));
	for(int i=1;i<=NumRows(A);i++)
		for(int j=1;j<=NumCols(A);j++) C[i][j]=(int)A[i][j];
	return C;
}


template <typename T> Matrix<float> Float(const Matrix<T>& A) {
	Matrix<float> C(NumRows(A),NumCols(A));
	for(int i=1;i<=NumRows(A);i++)
		for(int j=1;j<=NumCols(A);j++) C[i][j]=(float)A[i][j];
	return C;
}


template <typename T> Matrix<double> Double(const Matrix<T>& A) {
	Matrix<double> C(NumRows(A),NumCols(A));
	for(int i=1;i<=NumRows(A);i++)
		for(int j=1;j<=NumCols(A);j++) C[i][j]=(double)A[i][j];
	return C;
}


template <typename T> Vector<T> Moment(const Matrix<T>& A, int order){
	Matrix<T> C=A;
	Vector<T> mean=Mean(C,1);
	for(int i=1;i<=C.m;i++){
		C[i]-=mean;
		C[i]=Pow(C[i],order);
	}
	Vector<T> mom=Mean(C,1);
	return mom;
}


template <typename T> Vector<T> Var(const Matrix<T>& A,int k=0){
	Matrix<T> C=A;
	Vector<T> mean=Mean(C,1);
	for(int i=1;i<=NumRows(C);i++){
		C[i]-=mean;
		C[i]=Pow(C[i],2);
	}
	if(k==1) return Sum(C,1)/T(NumRows(C));
	else return Sum(C,1)/T(NumRows(C)-1);
}


template <typename T> inline Vector<T> Std(const Matrix<T>& A){
	return Pow(Var(A),0.5);
}


template <typename T> Vector<T> Skewness(const Matrix<T>& A){
	Vector<T> mom3=Moment(A,3);
	Vector<T> std3=Pow(Moment(A,2),1.5);
	return ArrayDivide(mom3,std3);
}


template <typename T> Vector<T> Kurtosis(const Matrix<T>& A){
	Vector<T> mom4=Moment(A,4);
	Vector<T> std4=Pow(Moment(A,2),2);
	return ArrayDivide(mom4,std4);
}


template <typename T> Matrix<T> Reshape(const Matrix<T>& A, int new_m, int new_n){
	assert(new_m*new_n == NumRows(A)*NumCols(A));
	Matrix<T> C(new_m,new_n);
	int i,j;
	int rowN=NumRows(A);
	for(j=1;j<=new_n;j++){
		for(i=1;i<=new_m;i++){
			int k=(j-1)*new_m+i-1;
			int colX=k/rowN+1;
			int rowX=k%rowN+1;
			C[i][j]=A[rowX][colX];
		}
	}
	return C;
}


template <typename T> double Norm(const Matrix<T>& A, int p=2) {
	int m=NumRows(A);
	int n=NumCols(A);
	int i,j;
	if(p==1){
		double maxsumabs=0;
		for(j=1;j<=n;j++){
			double sumabs=0;
			for(i=1;i<=m;i++) sumabs+=Abs(A[i][j]);
			if(maxsumabs<sumabs) maxsumabs=sumabs;
		}
		return maxsumabs;
	}
	else if(p==2){
		if(m==n){ return Matrix<T>::normnxn(A); }
		int l=local_min(m,n);
		Matrix<T> B(l,l);
		Matrix<T> trA=Transpose(A);
		if(m<n)
			B=A*trA;
		else
			B=trA*A;
		return sqrt(Matrix<T>::normnxn(B));
	}
	else if(p==mtl_numeric_limits<int>::max()){
		double maxsumabs=0;
		for(i=1;i<=m;i++) {
			double sumabs=0;
			for(j=1;j<=n;j++) sumabs+=Abs(A[i][j]);
			if(maxsumabs<sumabs) maxsumabs=sumabs;
		}
		return maxsumabs;
	}
	else{
		cerr << "ERROR: Invalid p-norm\n";
		assert(0);
		return 0;
	}
}


template <typename T> double NormF(const Matrix<T>& A) {
	return sqrt(real(Sum(Diag(Transpose(A)*A))));
}


template <typename T> double Matrix<T>::normnxn(const Matrix<T>& A) {
	int m=A.m; int n=A.n;
	Matrix<T> U(m,m), S(m,n), V(n,n);
	int a1=Svd(A,U,S,V);
	if(a1==0) {
		cerr << "Svd failed\n";
		assert(0);
		return 0;
	}
	double maxW=0;
	for(int i=1;i<=local