1 求解线性拟合模型如下的a1,a2,a3
已知 2a1+3a2+a3=4
4a1+2a2+3a1=17
7a1+a2-a3=-1
2 计算原理
A=xlsread('C:\Users\johon.ye\Desktop\x.xlsx','Sheet1');
B=xlsread('C:\Users\johon.ye\Desktop\y.xlsx','Sheet1');
M2=A.';
M3=M2*A;
M4=inv(M3);
M5=M4*M2;
M6=M5*B;
方法1 采用Opencv3.0中的Mat数据结构直接进行运算
程序代码如下：
#include "iostream"
#include "opencv2/core/core.hpp"
using namespace std;
using namespace cv;
int main(int argc, char *argv[])
{
Mat A = Mat::ones(3, 3, CV_32FC1);//输入变量X
Mat B = Mat::ones(3, 1, CV_32FC1);//输出变量Y
Mat X ;//输出方程求解结果
Mat m3;
Mat m4;
Mat m5;
Mat m6;
Mat m7;
A.at<float>(0, 0) = 2;
A.at<float>(0, 1) = 3;
A.at<float>(0, 2) = 1;
A.at<float>(1, 0) = 4;
A.at<float>(1, 1) = 2;
A.at<float>(1, 2) = 3;
A.at<float>(2, 0) = 7;
A.at<float>(2, 1) = 1;
A.at<float>(2, 2) = -1;
//A.at<float>(3, 0) = 3;
//A.at<float>(3, 1) = 2;
//A.at<float>(3, 2) = 4;
B.at<float>(0, 0) = 4;
B.at<float>(1, 0) = 17;
B.at<float>(2, 0) = 1;
//B.at<float>(3, 0) = 21;
X = ((A.t())*A).inv()*A.t()*B;
m3 = A.t();
m4 =( A.t())*A;
m5 = m4.inv();
m6 = m5*m3;
m7 = m6*B;
cout << "A=\n" << A << endl << endl;
cout << "B=\n" << B << endl << endl;
cout << "X=\n" <<X << endl << endl;
cout << "m3=\n" << m3 << endl << endl;
cout << "m4=\n" << m4 << endl << endl;
cout << "m5=\n" << m5 << endl << endl;
cout << "m6=\n" << m6<< endl << endl;
cout << "m7=\n" << m7<< endl << endl;
system("pause");
}
可以看出与实际原始结果有一定的偏差，主要原因是矩阵求逆运算过程引起的。增加一组
3a1+2a2+4a1=21 
增加一组1a1+2a2+a3=-19
通过增加X,y的取值查看是否可以优化计算发现输入的参数越多由于矩阵逆运算引起的偏差越大。
2 利用Matlab进行矩阵的运算
2a1+3a2+a3=4
4a1+2a2+3a1=17
7a1+a2-a3=-1程序如下
A=xlsread('C:\Users\johon.ye\Desktop\x.xlsx','Sheet1');
B=xlsread('C:\Users\johon.ye\Desktop\y.xlsx','Sheet1');
M2=A.';
M3=M2*A;
M4=inv(M3);
M5=M4*M2;
M6=M5*B;计算结果为
2a1+3a2+a3=4
4a1+2a2+3a1=17
7a1+a2-a3=-13a1+2a2+4a3=21
1a1+2a2+4a3=19
可见参数增加后，第二次计算结果优于第一次计算结果，因为系统一定不保留小数点点0，确定该结果正确，也可以看出Matlab有关计算运算的精度算法其整体优于Opencv中的Mat
方法3
使用c#依据矩阵运算的规则，编写矩阵类以及相应的运算函数。矩阵类的参考大牛代码，忘记出处了
//矩阵数据结构
//二维矩阵
class _Matrix
{
public int m;
public int n;
public float[] arr;
//初始化
public _Matrix()
{
m = 0;
n = 0;
}
public _Matrix(int mm, int nn)
{
m = mm;
n = nn;
}
//设置m
public void set_mn(int mm, int nn)
{
m = mm;
n = nn;
}
//设置m
public void set_m(int mm)
{
m = mm;
}
//设置n
public void set_n(int nn)
{
n = nn;
}
//初始化
public void init_matrix()
{
arr = new float[m * n];
}
//释放
public void free_matrix()
{
//delete [] arr;
}
//读取i,j坐标的数据
//失败返回-31415,成功返回值
public float read(int i, int j)
{
if (i >= m || j >= n)
{
return -31415;
}
//return *(arr + i * n + j);
return arr[i * n + j];
}
//写入i,j坐标的数据
//失败返回-1,成功返回1
public int write(int i, int j, float val)
{
if (i >= m || j >= n)
{
return -1;
}
arr[i * n + j] = val;
return 1;
}
}; //二维运算类
class _Matrix_Calc
{
//初始化
public _Matrix_Calc()
{
}
//C = A + B
//成功返回1,失败返回-1
public int add(ref _Matrix A, ref _Matrix B, ref _Matrix C)
{
int i = 0;
int j = 0;
//判断是否可以运算
if (A.m != B.m || A.n != B.n ||
A.m != C.m || A.n != C.n)
{
return -1;
}
//运算
for (i = 0; i < C.m; i++)
{
for (j = 0; j < C.n; j++)
{
C.write(i, j, A.read(i, j) + B.read(i, j));
}
}
return 1;
}
//C = A - B
//成功返回1,失败返回-1
public int subtract(ref _Matrix A, ref _Matrix B, ref _Matrix C)
{
int i = 0;
int j = 0;
//判断是否可以运算
if (A.m != B.m || A.n != B.n ||
A.m != C.m || A.n != C.n)
{
return -1;
}
//运算
for (i = 0; i < C.m; i++)
{
for (j = 0; j < C.n; j++)
{
C.write(i, j, A.read(i, j) - B.read(i, j));
}
}
return 1;
}
//C = A * B
//成功返回1,失败返回-1
public int multiply(ref _Matrix A, ref _Matrix B, ref _Matrix C)
{
int i = 0;
int j = 0;
int k = 0;
float temp = 0;
//判断是否可以运算
if (A.n != B.m)
{
return -1;
}
//运算
for (i = 0; i < C.m; i++)
{
for (j = 0; j < C.n; j++)
{
temp = 0;
for (k = 0; k < A.n; k++)
{
temp += A.read(i, k) * B.read(k, j);
}
C.write(i, j, temp);
}
}
return 1;
}
//行列式的值,只能计算2 * 2,3 * 3
//失败返回-31415,成功返回值
public float det(ref _Matrix A)
{
float value = 0;
//判断是否可以运算
if (A.m != A.n || (A.m != 2 && A.m != 3))
{
return -31415;
}
//运算
if (A.m == 2)
{
value = A.read(0, 0) * A.read(1, 1) - A.read(0, 1) * A.read(1, 0);
}
else
{
value = A.read(0, 0) * A.read(1, 1) * A.read(2, 2) +
A.read(0, 1) * A.read(1, 2) * A.read(2, 0) +
A.read(0, 2) * A.read(1, 0) * A.read(2, 1) -
A.read(0, 0) * A.read(1, 2) * A.read(2, 1) -
A.read(0, 1) * A.read(1, 0) * A.read(2, 2) -
A.read(0, 2) * A.read(1, 1) * A.read(2, 0);
}
return value;
}
//求转置矩阵,B = AT
//成功返回1,失败返回-1
public int transpos(ref _Matrix A, ref _Matrix B)
{
int i = 0;
int j = 0;
//判断是否可以运算
if (A.m != B.n || A.n != B.m)
{
return -1;
}
//运算
for (i = 0; i < B.m; i++)
{
for (j = 0; j < B.n; j++)
{
B.write(i, j, A.read(j, i));
}
}
return 1;
}
//求逆矩阵,B = A^(-1)
//成功返回1,失败返回-1
public int inverse(ref _Matrix A, ref _Matrix B)
{
int i = 0;
int j = 0;
int k = 0;
_Matrix m = new _Matrix(A.m, 2 * A.m);
float temp = 0;
float b = 0;
//判断是否可以运算
if (A.m != A.n || B.m != B.n || A.m != B.m)
{
return -1;
}
/*
//如果是2维或者3维求行列式判断是否可逆
if (A.m == 2 || A.m == 3)
{
if (det(A) == 0)
{
return -1;
}
}
*/
//增广矩阵m = A | B初始化
m.init_matrix();
for (i = 0; i < m.m; i++)
{
for (j = 0; j < m.n; j++)
{
if (j <= A.n - 1)
{
m.write(i, j, A.read(i, j));
}
else
{
if (i == j - A.n)
{
m.write(i, j, 1);
}
else
{
m.write(i, j, 0);
}
}
}
}
//高斯消元
//变换下三角
for (k = 0; k < m.m - 1; k++)
{
//如果坐标为k,k的数为0,则行变换
if (m.read(k, k) == 0)
{
for (i = k + 1; i < m.m; i++)
{
if (m.read(i, k) != 0)
{
break;
}
}
if (i >= m.m)
{
return -1;
}
else
{
//交换行
for (j = 0; j < m.n; j++)
{
temp = m.read(k, j);
m.write(k, j, m.read(k + 1, j));
m.write(k + 1, j, temp);
}
}
}
//消元
for (i = k + 1; i < m.m; i++)
{
//获得倍数
b = m.read(i, k) / m.read(k, k);
//行变换
for (j = 0; j < m.n; j++)
{
temp = m.read(i, j) - b * m.read(k, j);
m.write(i, j, temp);
}
}
}
//变换上三角
for (k = m.m - 1; k > 0; k--)
{
//如果坐标为k,k的数为0,则行变换
if (m.read(k, k) == 0)
{
for (i = k + 1; i < m.m; i++)
{
if (m.read(i, k) != 0)
{
break;
}
}
if (i >= m.m)
{
return -1;
}
else
{
//交换行
for (j = 0; j < m.n; j++)
{
temp = m.read(k, j);
m.write(k, j, m.read(k + 1, j));
m.write(k + 1, j, temp);
}
}
}
//消元
for (i = k - 1; i >= 0; i--)
{
//获得倍数
b = m.read(i, k) / m.read(k, k);
//行变换
for (j = 0; j < m.n; j++)
{
temp = m.read(i, j) - b * m.read(k, j);
m.write(i, j, temp);
}
}
}
//将左边方阵化为单位矩阵
for (i = 0; i < m.m; i++)
{
if (m.read(i, i) != 1)
{
//获得倍数
b = 1 / m.read(i, i);
//行变换
for (j = 0; j < m.n; j++)
{
temp = m.read(i, j) * b;
m.write(i, j, temp);
}
}
}
//求得逆矩阵
for (i = 0; i < B.m; i++)
{
for (j = 0; j < B.m; j++)
{
B.write(i, j, m.read(i, j + m.m));
}
}
//释放增广矩阵
m.free_matrix();
return 1;
}
};运算程序
static void Main(string[] args)
{
int i = 0;
int j = 0;
int k = 0;
_Matrix_Calc m_c = new _Matrix_Calc();
_Matrix m1 = new _Matrix(4, 3);
_Matrix m2 = new _Matrix(4, 1);
_Matrix m3 = new _Matrix(3, 4);//m1的转置
_Matrix m4=new _Matrix(3,3);//m1的转置与m1的乘积
_Matrix m5 = new _Matrix(3, 3);//m4的逆矩阵
_Matrix m6 = new _Matrix(3, 4);//m5与m3的乘积
_Matrix m7 = new _Matrix(3, 1);//m6与m2的乘积
//初始化内存
m1.init_matrix();
m2.init_matrix();
m3.init_matrix();
m4.init_matrix();
m5.init_matrix();
m6.init_matrix();
m7.init_matrix();
//初始化数据
k = 1;
float[,] temp1 = { { 2, 3, 1 }, { 4, 2, 3 }, { 7, 1, -1 }, { 3, 2, 4 } };
for (i = 0; i < m1.m; i++)
{
for (j = 0; j < m1.n; j++)
{
m1.write(i, j, temp1[i,j]);
}
}
Console.WriteLine("m1");
DisplayMatrix(m1);
// [4;17;1;21];
float[,] temp2 = new float[4, 1] {{4},{17},{1},{21} };
for (i = 0; i < m2.m; i++)
{
for (j = 0; j < m2.n; j++)
{
m2.write(i, j, temp2[i,j]);
}
}
Console.WriteLine("m2");
DisplayMatrix(m2);
m_c.transpos(ref m1, ref m3);
Console.WriteLine("m3");
DisplayMatrix(m3);
m_c.multiply(ref m3, ref m1, ref m4);
Console.WriteLine("m4");
DisplayMatrix(m4);
m_c.inverse(ref m4, ref m5);
Console.WriteLine("m5");
DisplayMatrix(m5);
m_c.multiply(ref m5, ref m3,ref m6);
Console.WriteLine("m6");
DisplayMatrix(m6);
m_c.multiply(ref m6, ref m2, ref m7);
Console.WriteLine("m7");
DisplayMatrix(m7);
}
static void DisplayMatrix(_Matrix DisMat )
{
int i;
int j;
for (i = 0; i < DisMat.m; i++)
{
for (j = 0; j < DisMat.n; j++)
{
string temp3= DisMat.read(i, j).ToString();
if (j <DisMat.n-1 )
{
Console.Write(temp3 + "\0");
}
else if (j == DisMat.n-1 )
{
Console.Write(temp3 + "\n");
}
}
}
//\n转意字符含义为换行
//
}
}
基本可以满足要求。
进一步测试发现存在问题
static void Main(string[] args)
{
int i = 0;
int j = 0;
int k = 0;
int Xrow = 8;
int Xcol=8;
int Yrow=8;
int Ycol = 1;
_Matrix_Calc m_c = new _Matrix_Calc();
_Matrix m1 = new _Matrix(Xrow, Xcol);
_Matrix m2 = new _Matrix(Yrow, Ycol);
_Matrix m3 = new _Matrix(Xcol, Xrow);//m1的转置
_Matrix m4 = new _Matrix(Xcol, Xcol);//m1的转置与m1的乘积
_Matrix m5 = new _Matrix(Xcol, Xcol);//m4的逆矩阵
_Matrix m6 = new _Matrix(Xcol, Xrow);//m5与m3的乘积
_Matrix m7 = new _Matrix(Xrow,Ycol);//m6与m2的乘积
//初始化内存
m1.init_matrix();
m2.init_matrix();
m3.init_matrix();
m4.init_matrix();
m5.init_matrix();
m6.init_matrix();
m7.init_matrix();
//初始化数据
k = 1;
//1.23 2.37 3.47 4.56 5.74 6.23 7.91 8.1
//1.24 2.34 3.46 4.55 5.73 6.21 7.9 8.11
//1.22 2.31 3.45 4.54 5.75 6.2 7.89 8.12
//1.21 2.3 3.44 4.53 5.76 6.22 7.88 8.13
//1.2 2.33 3.42 4.52 5.77 6.19 7.87 8.14
//1.27 2.305 3.43 4.51 5.78 6.18 7.86 8.15
//1.29 2.295 3.41 4.5 5.79 6.17 7.85 8.16
//1.19 2.345 3.4 4.49 5.72 6.16 7.93 8.17
float[,] temp1 = { {1.23f, 2.37f, 3.47f ,4.56f ,5.74f ,6.23f ,7.91f ,8.1f },
{ 1.24f ,2.34f, 3.46f, 4.55f, 5.73f, 6.21f, 7.9f, 8.11f },
{1.22f, 2.31f, 3.45f, 4.54f, 5.75f, 6.2f,7.89f, 8.12f},
{1.21f, 2.3f, 3.44f, 4.53f, 5.76f, 6.22f ,7.88f ,8.13f},
{ 1.2f, 2.33f, 3.42f, 4.52f, 5.77f, 6.19f, 7.87f ,8.14f },
{1.27f ,2.305f ,3.43f, 4.51f, 5.78f, 6.18f, 7.86f, 8.15f},
{1.29f ,2.295f, 3.41f, 4.5f, 5.79f, 6.17f, 7.85f, 8.16f},
{1.19f ,2.345f ,3.4f ,4.49f, 5.72f, 6.16f, 7.93f, 8.17f}
};
for (i = 0; i < m1.m; i++)
{
for (j = 0; j < m1.n; j++)
{
m1.write(i, j, temp1[i, j]);
}
}
Console.WriteLine("m1");
DisplayMatrix(m1);
// [4;17;1;21];
//43.895719
//43.895719
//43.815469
//43.750142
//43.737584
//43.711911
//43.759984
//43.739557
//43.6699
float[,] temp2 = { { 43.895719f }, { 43.815469f }, { 43.750142f }, { 43.737584f }, { 43.711911f }, { 43.759984f }, { 43.739557f }, { 43.6699f } };
for (i = 0; i < m2.m; i++)
{
for (j = 0; j < m2.n; j++)
{
m2.write(i, j, temp2[i, j]);
}
}
Console.WriteLine("m2");
DisplayMatrix(m2);
m_c.transpos(ref m1, ref m3);
Console.WriteLine("m3");
DisplayMatrix(m3);
m_c.multiply(ref m3, ref m1, ref m4);
Console.WriteLine("m4");
DisplayMatrix(m4);
m_c.inverse(ref m4, ref m5);
Console.WriteLine("m5");
DisplayMatrix(m5);
m_c.multiply(ref m5, ref m3, ref m6);
Console.WriteLine("m6");
DisplayMatrix(m6);
m_c.multiply(ref m6, ref m2, ref m7);
Console.WriteLine("m7");
DisplayMatrix(m7);
}
而通过matlab进行矩阵运算的结果与实际的结果相符合
同样的矩阵对于Opencv进行测试
扫一扫
4160
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
