计算数学
附录：矩阵及线性代数类设计方法
大部分高级语言并不是直接针对数值算法设计。比如C语言中，实现矩阵运算通常通过数组， 而二维动态数组实际通过指针实现。用户需要非常小心的分配与释放内存，否则内存泄漏造成程序异常。 这里以C++和C#为例，设计矩阵与线性代数类，通过类的设计，可以在析构函数中自动释放内存。 还可以通过索引器，使得C++或C#类似于MATLAB一样，实现矩阵运算。在java中不支持运算符重载， 因此java矩阵类，可以把各功能定义成函数。在Python等语言中，已经有第三方数值代数库支持。
C++矩阵向量类与线性代数
关于C++的语法书，请参考Stoustrup编写的《the C++ programming language》等专著或教材。头文件如下
/*------------------------------------------------------
Versions and Changes :
v1.0----2015-7-18
矩阵向量类与数值代数函数库
c++ matrix and vector library
主要包括
1. 简单的矩阵向量加减乘法运算，坐标旋转矩阵
2. 常用的矩阵分解 LDL, QR(MGS正交化，householder变换)
3. 最小二乘问题求解
4. 一般线性方程求解(正交化方法)
5. 矩阵求逆
6. 拉格朗日插值
Notes1:
1. 矩阵类的模板设计问题
矩阵向量类可以设计为模板，这样运算时候可以是的对各种数据类型具有通用性。
但这里仅把该类设计为double型，因为在数值运算中基本不可能使用单精度。
另外矩阵分解更不可能针对整型。实际运算中需要用整型矩阵运算的情况较少，
如果出现可以使用C++标准库vector构造，方法见note2。
2. 用c++ STL标准库中vector构造矩阵方法
构造方法(int ,double ,float)，以double为例
(1) vector > A(M,vector(N))
(2) vector > A(M)
for(int i=0;i(N)
(3) vector > A(N)
for(int i=0;i
Copyrigt(C) : Shanghai Astronomical Observatory, CAS
(All rights reserved) 2015
-------------------------------------------------------*/
#ifndef cmatlib_h
#define cmatlib_h
#include
#include
#include //format control
using namespace std;
class VEC
/*------------------------------------------------------
Versions and Changes :
v1.0----2015-7-18
vector class
下标索引算符为()
--------------------------------------------------------
Public Paras:
int size1 --------- rows of matrix
int size2 --------- cols of matri
Funcitons :
VEC (M) ---- creat a vector with lenght of M
operator() --- access the component
--------------------------------------------------------
Author : Song Yezhi
Copyrigt(C) : Shanghai Astronomical Observatory, CAS
(All rights reserved) 2015
-------------------------------------------------------*/
{
private:
double *x ;
int M1;
public:
// Constructors
VEC (int M)
//creat a vector
{
M1=M;
x= new double [M] ;
for(int i=0;i< M;i++)
x[i]=0.0;
}
// Destructor
~ VEC ()
{ delete [] x; }
double& operator() (int i)
// Component access
{ return x[i]; }
double operator () (int i) const { return x[i]; };
int size() const { return M1; };
void output(int width,int precision);
//output to the screen
void setv(double value);
//set the elements of the vector to a double value
//+,-,*,/ for the same index of two vectors
//which means that the size of the vectors must be the same
friend VEC operator + (const VEC& V1, const VEC& V2); //V1+V2
friend VEC operator - (const VEC& V1, const VEC& V2); //V1-V2
friend VEC operator * (const VEC& V1, const VEC& V2); //V1.*V2, for the matlab notation
friend VEC operator / (const VEC& V1, const VEC& V2); //V1./V2, for the matlab notation
//V1./V2
//+,-,*,/ for the vector and a scala
// the operation works on every elements of the vector
friend VEC operator + (const VEC& V1, double a); //V1+a
friend VEC operator - (const VEC& V1, double a); //V1-a
friend VEC operator * (const VEC& V1, double a); //V1*a
friend VEC operator / (const VEC& V1, double a); //V1/a
} ;
class MAT
/*------------------------------------------------------
Versions and Changes :
v1.0----2015-7-18
matrix class
在数组基础上构造,下标算符重载()
--------------------------------------------------------
Public Paras:
int size1 --------- rows of matrix
int size2 --------- cols of matrix
Functions :
MAT(M,N) ---- creat a vector with the size of (M,N)
operator() --- caccess the component
--------------------------------------------------------
Author : Song Yezhi
Copyrigt(C) : Shanghai Astronomical Observatory, CAS
(All rights reserved) 2015
-------------------------------------------------------*/
{
private:
double ** A;
int M1,N1;
public:
// Constructors
MAT(int M,int N)
{
M1=M;
N1=N;
int i,j;
A= new double *[M];
for(i=0;i< M;i++)
A[i]=new double [N] ;
// set the initial value to zero
for (i=0;i< M;i++)
for(j=0;j< N;j++)
A[i][j]=0.0;
}
// Destructor
~MAT()
{
int i;
for(i=0;i< M1;i++)
delete[] A[i];
delete[] A;
}
//Component access
double operator() (int i,int j) const { return A[i][j]; };
double & operator() (int i,int j) { return A[i][j]; };
int size1() const { return M1; };
int size2() const { return N1; };
void output(int width,int precision);
//output to the screen
void setv(double value);
//set the elements of the matrix to a double value
//+,-,*,/ for the same index of two matrixs
//which means that the size of the matrixs must be the same
friend MAT operator + (const MAT& A1, const MAT& A2); //A1+A2
friend MAT operator - (const MAT& A1, const MAT& A2); //A1-A2
friend MAT operator * (const MAT& A1, const MAT& A2); //A1.*A2
friend MAT operator / (const MAT& A1, const MAT& A2); //A1./A2
//+,-,*,/ for the matrix and a scala
// the operation works on every elements of the matrix
friend MAT operator + (const MAT& A1, double a); //A1+a
friend MAT operator - (const MAT& A1, double a); //A1-a
friend MAT operator * (const MAT& A1, double a); //A1*a
friend MAT operator / (const MAT& A1, double a); //A1/a
};
//--------------基本向量运算
void veccopy(VEC &b,const VEC &a);
//按值复制一个向量
double vecdot(const VEC &v1,const VEC & v2);
//向量内积
double norm(const VEC &v);
//向量2范数
//-----------------基本矩阵运算
void matcopy(MAT &B,const MAT &A);
//矩阵复制
void transpose(MAT &AT,const MAT &A);
//矩阵转置
void matmul(VEC &b,const MAT &A,const VEC &x);
// 矩阵乘以向量
void matmul(MAT &C,const MAT &A,const MAT &B);
//矩阵乘以矩阵
//------------------矩阵分解与求逆
void LDL(const MAT &A,MAT &L,VEC &D);
// LDL 分解
void MGS(const MAT &A, MAT &Q, MAT &R);
//修正的Gram-Schmidt正交化方法
void householder(const MAT &A,MAT &Q,MAT &R);
//hoseholder正交变换
void invlowtri(const MAT &R,MAT &S);
//inverse of lower triangula matrix
void invuptri(const MAT &U,MAT &R);
//inverse of upper triangula matrix
void invmat(const MAT &A,MAT &invA);
//对称正定矩阵逆矩阵
void inv(const MAT &A,MAT &iA);
//一般矩阵的逆矩阵 采用MGS分解方法
//--------------------线性代数方程
void LS_LDL(VEC &x,const MAT &A,const VEC &b);
//基于不开平方的cholesky分解计算最小二乘问题， A为对称正定矩阵
//void LS_hous(VEC &x,const MAT &A,const VEC &b);
//基于Householder变换求解最新二乘问题或适定问题
void LS_MGS( VEC &x,const MAT &A,const VEC &b);
//通过修正的Gram-Schimdt正交化求解最小二乘问题或适定问题
void LS_hous( VEC &x,const MAT &A,const VEC &b);
//least square solution or general linear equation
void uptri(const MAT &A,const VEC &b,VEC &x);
//上三角矩阵方程计算
void lowtri(const MAT &A, const VEC &b, VEC &x);
//下三角矩阵方程计算
void downtri(const MAT &A, const VEC &b, VEC &x);
//------------------插值算法
void lagrange(const VEC &x,const VEC &y,const VEC &xx,VEC &yy);
//lagrange interp
void rot_mat(MAT &Rmat,double angle,char axID);
//rotation matrix
double robustweight(double v,double sigma);
//used for robust estimation
#endif
C++代码如下：
/*------------------------------------------------------
Versions and Changes :
v1.0----2015-7-18
c++ matrix library
the comment is in cmatlib.h
v2.0---2015-11-05
the second editon of cmatlib
--------------------------------------------------------
Notes:
1. 头文件为 matlib.h
2. 调用该函数库需要加入头文件
#include "matlib.h"
3. 编译方式为 $(cc) main.cpp matlib.cpp
--------------------------------------------------------
Author : Song Yezhi
Copyrigt(C) : Shanghai Astronomical Observatory, CAS
(All rights reserved) 2015
-------------------------------------------------------*/
//#include
//#include
//#include //
#include "cmatlib.h"
//
//using namespace std;
void VEC::output(int width,int precision)
/*------------------------------------------------------
Versions and Changes :
v1.0---- 2015-09-11
output the vector to the screen
--------------------------------------------------------
Input Parameters :
width --- the width of the element value
precision-- presicion of the element value
Output Parameters :
Notes :
--------------------------------------------------------
Author : Song Yezhi
Copyrigt(C) : Shanghai Astronomical Observatory, CAS
(All rights reserved) 2015
-------------------------------------------------------*/
{
int i;
for(i=0;i< M1;i++)
cout<
Copyrigt(C) : Shanghai Astronomical Observatory, CAS
(All rights reserved) 2015
-------------------------------------------------------*/
{
int i;
for(i=0;i< M1;i++)
x[i]=value;
}
void MAT::output(int width,int precision)
/*------------------------------------------------------
Versions and Changes :
v1.0---- 2015-09-11
output the matrix to the screen
in many case this function is used for debug
--------------------------------------------------------
Input Parameters :
width --- the width of the element value
precision-- presicion of the element value
Output Parameters :
Notes :
--------------------------------------------------------
Author : Song Yezhi
Copyrigt(C) : Shanghai Astronomical Observatory, CAS
(All rights reserved) 2015
-------------------------------------------------------*/
{
int i,j;
for(i=0;i< M1;i++)
{
for(j=0;j< N1;j++) cout<
Copyrigt(C) : Shanghai Astronomical Observatory, CAS
(All rights reserved) 2015
-------------------------------------------------------*/
{
int i,j;
for(i=0;i< M1;i++)
{
for(j=0;j< N1;j++)
A[i][j]=value;
}
}
VEC operator + (const VEC& V1, const VEC& V2)
/*------------------------------------------------------
Versions and Changes :
v1.0---- 2015-11-05
addition of the the two vectors
V1+V2
--------------------------------------------------------
Input Parameters :
V1
V2
Output Parameters :
V3=
Notes :
--------------------------------------------------------
Author : Song Yezhi
Copyrigt(C) : Shanghai Astronomical Observatory, CAS
(All rights reserved) 2015
-------------------------------------------------------*/
{
int N=V1.size();
int i;
VEC VT(N);
for(i=0;i< N;i++)
VT(i)=V1(i)+V2(i);
return VT;
}
VEC operator - (const VEC& V1, const VEC& V2)
/*------------------------------------------------------
Versions and Changes :
v1.0---- 2015-11-05
vector subtraction V1-V2
--------------------------------------------------------
Input Parameters :
Output Parameters :
Notes :
--------------------------------------------------------
Author : Song Yezhi
Copyrigt(C) : Shanghai Astronomical Observatory, CAS
(All rights reserved) 2015
-------------------------------------------------------*/
{
int N=V1.size();
int i;
VEC VT(N);
for(i=0;i< N;i++)
VT(i)=V1(i)-V2(i);
return VT;
}
VEC operator * (const VEC& V1, const VEC& V2)
/*------------------------------------------------------
Versions and Changes :
v1.0---- 2015-11-05
V1.*V2
--------------------------------------------------------
Input Parameters :
Output Parameters :
Notes :
--------------------------------------------------------
Author : Song Yezhi
Copyrigt(C) : Shanghai Astronomical Observatory, CAS
(All rights reserved) 2015
-------------------------------------------------------*/
{
int N=V1.size();
int i;
VEC VT(N);
for(i=0;i< N;i++)
VT(i)=V1(i)*V2(i);
return VT;
}
VEC operator / (const VEC& V1, const VEC& V2)
/*------------------------------------------------------
Versions and Changes :
v1.0---- 2015-11-05
V1./V2
--------------------------------------------------------
Input Parameters :
Output Parameters :
Notes :
--------------------------------------------------------
Author : Song Yezhi
Copyrigt(C) : Shanghai Astronomical Observatory, CAS
(All rights reserved) 2015
-------------------------------------------------------*/
{
int N=V1.size();
int i;
VEC VT(N);
for(i=0;i< N;i++)
VT(i)=V1(i)/V2(i);
return VT;
}
VEC operator + (const VEC& V1,double a)
/*------------------------------------------------------
Versions and Changes :
v1.0---- 2015-11-05
V1+a
--------------------------------------------------------
Input Parameters :
Output Parameters :
Notes :
--------------------------------------------------------
Author : Song Yezhi
Copyrigt(C) : Shanghai Astronomical Observatory, CAS
(All rights reserved) 2015
-------------------------------------------------------*/
{
int N=V1.size();
int i;
VEC VT(N);
for(i=0;i< N;i++)
VT(i)=V1(i)+a;
return VT;
}
VEC operator - (const VEC& V1,double a)
/*------------------------------------------------------
Versions and Changes :
v1.0---- 2015-11-05
V1-a
--------------------------------------------------------
Input Parameters :
Output Parameters :
Notes :
--------------------------------------------------------
Author : Song Yezhi
Copyrigt(C) : Shanghai Astronomical Observatory, CAS
(All rights reserved) 2015
-------------------------------------------------------*/
{
int N=V1.size();
int i;
VEC VT(N);
for(i=0;i< N;i++)
VT(i)=V1(i)-a;
return VT;
}
VEC operator * (const VEC& V1,double a)
/*------------------------------------------------------
Versions and Changes :
v1.0---- 2015-11-05
V1*a
--------------------------------------------------------
Input Parameters :
Output Parameters :
Notes :
--------------------------------------------------------
Author : Song Yezhi
Copyrigt(C) : Shanghai Astronomical Observatory, CAS
(All rights reserved) 2015
-------------------------------------------------------*/
{
int N=V1.size();
int i;
VEC VT(N);
for(i=0;i< N;i++)
VT(i)=V1(i)*a;
return VT;
}
VEC operator / (const VEC& V1,double a)
/*------------------------------------------------------
Versions and Changes :
v1.0---- 2015-11-05
V1/a
--------------------------------------------------------
Input Parameters :
Output Parameters :
Notes :
--------------------------------------------------------
Author : Song Yezhi
Copyrigt(C) : Shanghai Astronomical Observatory, CAS
(All rights reserved) 2015
-------------------------------------------------------*/
{
int N=V1.size();
int i;
VEC VT(N);
for(i=0;i< N;i++)
VT(i)=V1(i)/a;
return VT;
}
void veccopy(VEC &b,const VEC &a)
/*------------------------------------------------------
Versions and Changes :
v1.0----2015-8-5
按值复制一个向量
--------------------------------------------------------
Input Parameters :
a
Output Parameters :
b
--------------------------------------------------------
Author : Song Yezhi
Copyrigt(C) : Shanghai Astronomical Observatory, CAS
(All rights reserved) 2015
-------------------------------------------------------*/
{
int N=a.size();
int M=b.size();
if(N!=M)
{cout<<"the dimention of two vectors must be the same"<
Copyrigt(C) : Shanghai Astronomical Observatory, CAS
(All rights reserved) 2015
------------------------------------------------------- */
{
int N=v1.size();
double sum=0;
for(int i=0;i< N;i++)
sum=sum+v1(i)*v2(i);
return sum;
}
double norm(const VEC &v)
/*------------------------------------------------------
Versions and Changes :
v1.0----2015-7-30
norm-2 of a vector
--------------------------------------------------------
Input Parameters :
v1--- vector
v2--- vector
Output Parameters :
a --- double
--------------------------------------------------------
Author : Song Yezhi
Copyrigt(C) : Shanghai Astronomical Observatory, CAS
(All rights reserved) 2015
------------------------------------------------------- */
{
int N=v.size();
double sum=0;
for(int i=0;i< N;i++)
sum=sum+v(i)*v(i);
sum=sqrt(sum);
return sum;
}
void matcopy(MAT &B,const MAT &A)
/* ------------------------------------------------------
Versions and Changes :
v1.0----2015-8-5
按值复制矩阵
MAT matcopy(const MAT &A)
v1.1--- 改为void类型函数
通过地址传递变了，因为MAT的()算符重载问题
--------------------------------------------------------
Author : Song Yezhi
Copyrigt(C) : Shanghai Astronomical Observatory, CAS
(All rights reserved) 2015
------------------------------------------------------- */
{
int M=A.size1();
int N=A.size2();
int i,j;
for(i=0;i< M;i++)
for(j=0;j< N;j++)
B(i,j)=A(i,j);
}
void transpose(MAT &AT,const MAT &A)
/* ------------------------------------------------------
Versions and Changes :
v1.0----2015-7-16
矩阵转置
--------------------------------------------------------
Author : Song Yezhi
Copyrigt(C) : Shanghai Astronomical Observatory, CAS
(All rights reserved) 2015
------------------------------------------------------- */
{
int M, N;
M = A.size1();
N = A.size2();
int P=AT.size1();
int Q=AT.size2();
if(N!=P||M!=Q)
{cout<<"the dimention of two vectors must be the same"<
Copyrigt(C) : Shanghai Astronomical Observatory, CAS
(All rights reserved) 2015
-------------------------------------------------------*/
{
int row=A1.size1();
int col=A1.size2();
MAT A3(row,col);
for(int i=0;i< row;i++)
for(int j=0;j< col;j++)
A3(i,j)=A1(i,j)+A2(i,j);
return A3;
}
MAT operator -(const MAT& A1,const MAT& A2)
/*------------------------------------------------------
Versions and Changes :
v1.0---- 2015-11-05
A1-A2
--------------------------------------------------------
Input Parameters :
Output Parameters :
Notes :
--------------------------------------------------------
Author : Song Yezhi
Copyrigt(C) : Shanghai Astronomical Observatory, CAS
(All rights reserved) 2015
-------------------------------------------------------*/
{
int row=A1.size1();
int col=A1.size2();
MAT A3(row,col);
for(int i=0;i< row;i++)
for(int j=0;j< col;j++)
A3(i,j)=A1(i,j)-A2(i,j);
return A3;
}
MAT operator *(const MAT& A1,const MAT& A2)
/*------------------------------------------------------
Versions and Changes :
v1.0---- 2015-11-05
A1.*A2
--------------------------------------------------------
Input Parameters :
Output Parameters :
Notes :
--------------------------------------------------------
Author : Song Yezhi
Copyrigt(C) : Shanghai Astronomical Observatory, CAS
(All rights reserved) 2015
-------------------------------------------------------*/
{
int row=A1.size1();
int col=A1.size2();
MAT A3(row,col);
for(int i=0;i< row;i++)
for(int j=0;j< col;j++)
A3(i,j)=A1(i,j)*A2(i,j);
return A3;
}
MAT operator /(const MAT& A1,const MAT& A2)
/*------------------------------------------------------
Versions and Changes :
v1.0---- 2015-11-05
A1./A2
--------------------------------------------------------
Input Parameters :
Output Parameters :
Notes :
--------------------------------------------------------
Author : Song Yezhi
Copyrigt(C) : Shanghai Astronomical Observatory, CAS
(All rights reserved) 2015
-------------------------------------------------------*/
{
int row=A1.size1();
int col=A1.size2();
MAT A3(row,col);
for(int i=0;i< row;i++)
for(int j=0;j< col;j++)
A3(i,j)=A1(i,j)/A2(i,j);
return A3;
}
MAT operator +(const MAT& A1,double a)
/*------------------------------------------------------
Versions and Changes :
v1.0---- 2015-11-05
A1+a
--------------------------------------------------------
Input Parameters :
Output Parameters :
Notes :
--------------------------------------------------------
Author : Song Yezhi
Copyrigt(C) : Shanghai Astronomical Observatory, CAS
(All rights reserved) 2015
-------------------------------------------------------*/
{
int row=A1.size1();
int col=A1.size2();
MAT A3(row,col);
for(int i=0;i< row;i++)
for(int j=0;j< col;j++)
A3(i,j)=A1(i,j)+a;
return A3;
}
MAT operator -(const MAT& A1,double a)
/*------------------------------------------------------
Versions and Changes :
v1.0---- 2015-11-05
A1-a
--------------------------------------------------------
Input Parameters :
Output Parameters :
Notes :
--------------------------------------------------------
Author : Song Yezhi
Copyrigt(C) : Shanghai Astronomical Observatory, CAS
(All rights reserved) 2015
-------------------------------------------------------*/
{
int row=A1.size1();
int col=A1.size2();
MAT A3(row,col);
for(int i=0;i< row;i++)
for(int j=0;j< col;j++)
A3(i,j)=A1(i,j)-a;
return A3;
}
MAT operator *(const MAT& A1,double a)
/*------------------------------------------------------
Versions and Changes :
v1.0---- 2015-11-05
A1*a
--------------------------------------------------------
Input Parameters :
Output Parameters :
Notes :
--------------------------------------------------------
Author : Song Yezhi
Copyrigt(C) : Shanghai Astronomical Observatory, CAS
(All rights reserved) 2015
-------------------------------------------------------*/
{
int row=A1.size1();
int col=A1.size2();
MAT A3(row,col);
for(int i=0;i< row;i++)
for(int j=0;j< col;j++)
A3(i,j)=A1(i,j)*a;
return A3;
}
MAT operator /(const MAT& A1,double a)
/*------------------------------------------------------
Versions and Changes :
v1.0---- 2015-11-05
A1/a
--------------------------------------------------------
Input Parameters :
Output Parameters :
Notes :
--------------------------------------------------------
Author : Song Yezhi
Copyrigt(C) : Shanghai Astronomical Observatory, CAS
(All rights reserved) 2015
-------------------------------------------------------*/
{
int row=A1.size1();
int col=A1.size2();
MAT A3(row,col);
for(int i=0;i< row;i++)
for(int j=0;j< col;j++)
A3(i,j)=A1(i,j)/a;
return A3;
}
void matmul(VEC &b,const MAT &A,const VEC &x)
/*-----------------------------------------------------
Versions and Changes :
v1.0----2015-7-17
b=Ax
--------------------------------------------------------
Input Parameters :
A(M,N)--- matrix
x(N)--- vector
Output Parameters :
b(M) --- vector
--------------------------------------------------------
Author : Song Yezhi
Copyrigt(C) : Shanghai Astronomical Observatory, CAS
(All rights reserved) 2015
-------------------------------------------------------*/
{
int M,N;
M=A.size1();
N=A.size2();
double sum;
int i,j;
for (i=0;i< M;i++)
{
sum=0.0;
for(j=0;j< N;j++)
sum=sum+A(i,j)*x(j);
b(i)=sum;
}
}
void matmul(MAT &C,const MAT &A,const MAT &B)
/* ------------------------------------------------------
Versions and Changes :
v1.0----2015-7-17
C=AB
--------------------------------------------------------
Input Parameters :
A(M,N)--- matrix
B(N,P)--- matrix
Output Parameters :
C(M,P) --- matrix
--------------------------------------------------------
Author : Song Yezhi
Copyrigt(C) : Shanghai Astronomical Observatory, CAS
(All rights reserved) 2015
-------------------------------------------------------*/
{
int M=A.size1();
int N=A.size2();
int P=C.size2();
int i,j,k;
double tmp1,tmp2;
double sum;
for(i=0;i< M;i++)
for(j=0;j< P;j++)
{
sum=0.0;
for(k=0;k< N;k++) sum=sum+A(i,k)*B(k,j);
C(i,j)=sum;
}
return;
}