矩阵操作

发布于 2021-01-04  150 次阅读


markdown latex测试

题目

编写二维实数矩阵(非方阵)的相关运算函数,并在主函数测试各个函数的功能,并输出计算结果。具体要求如下:
1. 实现两个矩阵的加减法: A+B, A-B (注意维度大小匹配)
2. 实现两个矩阵的乘法: A*B (注意维度大小匹配)
3. 求给定矩阵的转置矩阵: A^T
4. 计算给定方阵的迹(主对角元素之和) \mathrm{Tr}(A)
5. 计算给定方阵的主对角元素之积
6. 计算给定矩阵的1-范数 (最大的行绝对值和)
7. 计算给定矩阵的 无穷-范数 (最大的列绝对值和)
8. 计算给定矩阵的 弗罗贝尼乌斯范数(Frobenius)范数 (定义为所有矩阵元素的平方和再开根号;提示:也可以利用矩阵A乘以其自身的转置矩阵A^T然后取所得方阵的迹进行计算)

一、实现方法

1.1 实现思路

题目要求为实现矩阵本身和矩阵与矩阵之间的一些数值计算和操作。总体思想为,使用两个二维数组a、b来记录题目所给的矩阵,并用row1、col1、row2、col2来记录两个矩阵的行数和列数。 c为记录输出矩阵的二维数组。采用循环的方式实现题目所要求的各个操作。

二、函数实现

2.1 两个矩阵的加减法

注:由于题目未指明上限,此处对于数组的定义采用

 #define N 100

对于一个大小为100*100的数组,对于题目中所需要的全局数组,其大小

tot = 2 * sizeof(a);

为163216Byte,即0.15Mb,对于现代计算机来说非常小,而且可以满足绝大部分实际计算需要。

\begin{pmatrix}
a_{1,1}&a_{1,2}&a_{1,3}\\
a_{2,1}&a_{2,2}&a_{2,3}\\
a_{3,1}&a_{3,2}&a_{3,3}\\
\end{pmatrix}+
\begin{pmatrix}
b_{1,1}&b_{1,2}&b_{1,3}\\
b_{2,1}&b_{2,2}&b_{2,3}\\
b_{3,1}&b_{3,2}&b_{3,3}\\
\end{pmatrix}=
\begin{pmatrix}
c_{1,1}&c_{1,2}&c_{1,3}\\
c_{2,1}&c_{2,2}&c_{2,3}\\
c_{3,1}&c_{3,2}&c_{3,3}\\
\end{pmatrix}

\begin{pmatrix}
a_{1,1}&a_{1,2}&a_{1,3}\\
a_{2,1}&a_{2,2}&a_{2,3}\\
a_{3,1}&a_{3,2}&a_{3,3}\\
\end{pmatrix}-
\begin{pmatrix}
b_{1,1}&b_{1,2}&b_{1,3}\\
b_{2,1}&b_{2,2}&b_{2,3}\\
b_{3,1}&b_{3,2}&b_{3,3}\\
\end{pmatrix}=
\begin{pmatrix}
c_{1,1}&c_{1,2}&c_{1,3}\\
c_{2,1}&c_{2,2}&c_{2,3}\\
c_{3,1}&c_{3,2}&c_{3,3}\\
\end{pmatrix}

两个函数的加减法,采用一般的循环进行相加和相减,存入c中之后进行输出,代码如下。
其中,将a[][N+1],b[][N+1]作为传递数组,省去数组复制的时间,同时在函数中保留原数组的行数、列数信息,用于后续的操作。

void p(double a[][N + 1], double b[][N + 1], int row1, int col1, int row2, int col2) //矩阵加法
{
    double c[N + 1][N + 1];
    memset(c, 0, sizeof(c));
    if (row1 != row2 || col1 != col2)
    {
        cout << "Can't do this!" << endl;
        cout << endl;
        return;
    }
    for (int i = 1; i <= row1; i++)
    {
        for (int j = 1; j <= col1; j++)
                c[i][j] = a[i][j] + b[i][j]; //答案存储矩阵
    }
    for (int i = 1; i <= row1; i++)//输出
    {
        for (int j = 1; j <= col2; j++)
            cout << c[i][j] << ' ';
        cout << endl;
    }
    cout << endl;
}
void d(double a[][N + 1], double b[][N + 1], int row1, int col1, int row2, int col2) //矩阵减法
{
    double c[N + 1][N + 1];
    memset(c, 0, sizeof(c));
    if (row1 != row2 || col1 != col2)
    {
        cout << "Can't do this!" << endl;
        cout << endl;
        return;
    }
    for (int i = 1; i <= row1; i++)
    {
        for (int j = 1; j <= col1; j++)
            c[i][j] = a[i][j] - b[i][j]; //答案存储矩阵
    }
    for (int i = 1; i <= row1; i++)//输出
    {
        for (int j = 1; j <= col2; j++)
            cout << c[i][j] << ' ';
        cout << endl;
    }
    cout << endl;
}

首先,为了确保两个矩阵是可以进行相加和相减的,进行了一个判断,当两个矩阵大小一致才进行后续的操作。

2.2 两个矩阵的乘法

\begin{pmatrix}
a_{1,1}&a_{1,2}&a_{1,3}\\
a_{2,1}&a_{2,2}&a_{2,3}\\
\end{pmatrix}*
\begin{pmatrix}
b_{1,1}&b_{1,2}\\
b_{2,1}&b_{2,2}\\
b_{3,1}&b_{3,2}\\
\end{pmatrix}=
\begin{pmatrix}
c_{1,1}&c_{1,2}\\
c_{2,1}&c_{2,2}\\
\end{pmatrix}

矩阵的点乘,根据定义,采用三个for循环进行嵌套实现,复杂度为O(n^3),对于我们定义的数组大小来说,可以满足一般需要,如果对于更大的矩阵,可以采用Strassen演算法优化,或是采用并行矩阵乘法。
同样的,需要判断矩阵a的列数和矩阵b的行数是否一致再进行后续操作。

 void dot(double a[][N + 1], double b[][N + 1], int row1, int col1, int row2, int col2)//矩阵乘法
 {
    double c[N + 1][N + 1];
    memset(c, 0, sizeof(c));
    if (row2 != col1)
    {
        cout << "Can't do this!" << endl;
        cout << endl;
        return;
    }
    for (int i = 1; i <= row1; i++)
    {
        for (int j = 1; j <= col2; j++)
            for (int k = 1; k <= col1; k++)
                c[i][j] += a[i][k] * b[k][j]; //答案存储矩阵
    }
    for (int i = 1; i <= row1; i++)//输出
    {
        for (int j = 1; j <= col2; j++)
            cout << c[i][j] << ' ';
        cout << endl;
    }
    cout << endl;
 }

2.3 矩阵的转置

\begin{pmatrix}
a_{1,1}&a_{1,2}&a_{1,3}\\
a_{2,1}&a_{2,2}&a_{2,3}\\
\end{pmatrix}^T=
\begin{pmatrix}
a_{1,1}&a_{2,1}\\
a_{1,2}&a_{2,2}\\
a_{1,3}&a_{2,3}\\
\end{pmatrix}

void T(double a[][N + 1], int row1, int col1) //矩阵转置
{
    double c[N + 1][N + 1];
    memset(c, 0, sizeof(c));
    for (int i = 1; i <= row1; i++)
    {
        for (int j = 1; j <= col1; j++)
            c[j][i] = a[i][j];
    }
    for (int i = 1; i <= col1; i++)//输出
    {
        for (int j = 1; j <= row1; j++)
            cout << c[i][j] << ' ';
        cout << endl;
    }
    cout << endl;
}

2.4 矩阵的迹

tr
\begin{pmatrix}
a_{1,1}&a_{1,2}&a_{1,3}\\
a_{2,1}&a_{2,2}&a_{2,3}\\
a_{3,1}&a_{3,2}&a_{3,3}\\
\end{pmatrix}=
a_{1,1} + a_{2,2} + a_{3,3}


矩阵的迹的实现使用一个循环直接求对角数之和即可,此处的判断需要为一个方阵才可以进行求迹操作。

void tr(double a[][N + 1], int row1, int col1) //矩阵的迹
{
    if (row1 != col1)
    {
        cout << "Can't do this!" << endl;
        cout << endl;
        return;
    }
    double sum = 0;
    for (int i = 1; i <= row1; i++)
        sum += a[i][i];
    cout << sum << endl;
    cout << endl;
}

2.5 矩阵的对角元素之积

将该操作即为mul,则有

mul
\begin{pmatrix}
a_{1,1}&a_{1,2}&a_{1,3}\\
a_{2,1}&a_{2,2}&a_{2,3}\\
a_{3,1}&a_{3,2}&a_{3,3}\\
\end{pmatrix}=
a_{1,1} * a_{2,2} * a_{3,3}


类似于求迹,将加法换成乘法即可。此处需要特别注意的是,不可以将存储答案的变量初始化为0,这里采用了将其初始化为第一个元素的方法。同样需要判断是否为方阵。

void mul(double a[][N + 1], int row1, int col1) //矩阵对角线积
{
    if (row1 != col1)
    {
        cout << "Can't do this!" << endl;
        cout << endl;
        return;
    }
    double sum = a[1][1];
    for (int i = 2; i <= row1; i++)
        sum *= a[i][i];
    cout << sum << endl;
    cout << endl;
}

2.6 矩阵的1-范数

ans = max(\sum\limits_{j = 1}^{rol} |a_{i,j}|) (i = 1, 2,3, ······,row)


使用自行构造的求最大值函数,以简化程序。

double max(double x, double y) { return x > y ? x : y; }

求1-范数采用两个循环,求出每一行的绝对值之和再取最大值即可。

void F1(double a[][N + 1], int row1, int col1) //矩阵的1-范数
{
    double sum = 0, ans = 0;
    for (int i = 1; i <= row1; i++)
    {
        for (int j = 1; j <= col1; j++)
            sum += abs(a[i][j]);
        ans = max(ans, sum);
        sum = 0;
    }
    cout << ans << endl;
    cout << endl;
}

2.7 矩阵的无穷-范数

ans = max(\sum\limits_{j = 1}^{row} |a_{j,i}|) (i = 1, 2,3, ······,rol)


和1-范数求法类似,将i,j简单对调即可。

void F2(double a[][N + 1], int row1, int col1)//矩阵的无穷-范数
{
    double sum = 0, ans = 0;
    for (int i = 1; i <= col1; i++)
    {
        for (int j = 1; j <= row1; j++)
            sum += abs(a[j][i]);
        ans = max(ans, sum);
        sum = 0;
    }
    cout << ans << endl;
    cout << endl;
}

2.8 矩阵的Frobenius范数

ans=\sqrt{\sum\limits_{i=1}^{row}\sum\limits_{j=1}^{col}|a_{i,j}|^2}


直接求出每个元素的平方和再使用cmath中的sqrt函数求开方即可。

void frobenius(double a[][N + 1], int row1, int col1)//矩阵的Frobenius 范数
{
    double sum = 0;
    for (int i = 1; i <= row1; i++)
        for (int j = 1; j <= col1; j++)
            sum += a[i][j] * a[i][j];
    cout << sqrt(sum) << endl;
    cout << endl;
}

下面给出完整程序,并附上样例输入输出。

#include <stdio.h>
#include <iostream>
#include <cmath>
#include <string.h>
using namespace std;
#define N 100
double max(double x, double y) { return x > y ? x : y; }
double a[N + 1][N + 1], b[N + 1][N + 1];
int max1, max2;
void dot(double a[][N + 1], double b[][N + 1], int row1, int col1, int row2, int col2)//矩阵乘法
{
    double c[N + 1][N + 1];
    memset(c, 0, sizeof(c));
    if (row2 != col1)
    {
        cout << "Can't do this!" << endl;
        cout << endl;
        return;
    }
    for (int i = 1; i <= row1; i++)
    {
        for (int j = 1; j <= col2; j++)
            for (int k = 1; k <= col1; k++)
                c[i][j] += a[i][k] * b[k][j]; //答案存储矩阵
    }
    for (int i = 1; i <= row1; i++)//输出
    {
        for (int j = 1; j <= col2; j++)
            cout << c[i][j] << ' ';
        cout << endl;
    }
    cout << endl;
}
void p(double a[][N + 1], double b[][N + 1], int row1, int col1, int row2, int col2) //矩阵加法
{
    double c[N + 1][N + 1];
    memset(c, 0, sizeof(c));
    if (row1 != row2 || col1 != col2)
    {
        cout << "Can't do this!" << endl;
        cout << endl;
        return;
    }
    for (int i = 1; i <= row1; i++)
    {
        for (int j = 1; j <= col1; j++)
                c[i][j] = a[i][j] + b[i][j]; //答案存储矩阵
    }
    for (int i = 1; i <= row1; i++)//输出
    {
        for (int j = 1; j <= col2; j++)
            cout << c[i][j] << ' ';
        cout << endl;
    }
    cout << endl;
}
void d(double a[][N + 1], double b[][N + 1], int row1, int col1, int row2, int col2) //矩阵减法
{
    double c[N + 1][N + 1];
    memset(c, 0, sizeof(c));
    if (row1 != row2 || col1 != col2)
    {
        cout << "Can't do this!" << endl;
        cout << endl;
        return;
    }
    for (int i = 1; i <= row1; i++)
    {
        for (int j = 1; j <= col1; j++)
            c[i][j] = a[i][j] - b[i][j]; //答案存储矩阵
    }
    for (int i = 1; i <= row1; i++)//输出
    {
        for (int j = 1; j <= col2; j++)
            cout << c[i][j] << ' ';
        cout << endl;
    }
    cout << endl;
}
void T(double a[][N + 1], int row1, int col1) //矩阵转置
{
    double c[N + 1][N + 1];
    memset(c, 0, sizeof(c));
    for (int i = 1; i <= row1; i++)
    {
        for (int j = 1; j <= col1; j++)
            c[j][i] = a[i][j];
    }
    for (int i = 1; i <= col1; i++)//输出
    {
        for (int j = 1; j <= row1; j++)
            cout << c[i][j] << ' ';
        cout << endl;
    }
    cout << endl;
}
void tr(double a[][N + 1], int row1, int col1) //矩阵的迹
{
    if (row1 != col1)
    {
        cout << "Can't do this!" << endl;
        cout << endl;
        return;
    }
    double sum = 0;
    for (int i = 1; i <= row1; i++)
        sum += a[i][i];
    cout << sum << endl;
    cout << endl;
}
void mul(double a[][N + 1], int row1, int col1) //矩阵对角线积
{
    if (row1 != col1)
    {
        cout << "Can't do this!" << endl;
        cout << endl;
        return;
    }
    double sum = a[1][1];
    for (int i = 2; i <= row1; i++)
        sum *= a[i][i];
    cout << sum << endl;
    cout << endl;
}
void F1(double a[][N + 1], int row1, int col1) //矩阵的1-范数
{
    double sum = 0, ans = 0;
    for (int i = 1; i <= row1; i++)
    {
        for (int j = 1; j <= col1; j++)
            sum += abs(a[i][j]);
        ans = max(ans, sum);
        sum = 0;
    }
    cout << ans << endl;
    cout << endl;
}
void F2(double a[][N + 1], int row1, int col1)//矩阵的无穷-范数
{
    double sum = 0, ans = 0;
    for (int i = 1; i <= col1; i++)
    {
        for (int j = 1; j <= row1; j++)
            sum += abs(a[j][i]);
        ans = max(ans, sum);
        sum = 0;
    }
    cout << ans << endl;
    cout << endl;
}
void frobenius(double a[][N + 1], int row1, int col1)//矩阵的Frobenius 范数
{
    double sum = 0;
    for (int i = 1; i <= row1; i++)
        for (int j = 1; j <= col1; j++)
            sum += a[i][j] * a[i][j];
    cout << sqrt(sum) << endl;
    cout << endl;
}

int main()
{
    int row1, col1, row2, col2;
    cin >> row1 >> col1;
    for (int i = 1; i <= row1; i++)
        for (int j = 1; j <= col1; j++)
            cin >> a[i][j];
    cin >> row2 >> col2;
    for (int i = 1; i <= row2; i++)
        for (int j = 1; j <= col2; j++)
            cin >> b[i][j];
    printf("矩阵乘法\n");
    dot(a, b, row1, col1, row2, col2);
    printf("矩阵加法\n");
    p(a, b, row1, col1, row2, col2);
    printf("矩阵减法\n");
    d(a, b, row1, col1, row2, col2);
    printf("矩阵转置\n");
    T(a, row1, col1);
    printf("矩阵的迹\n");
    tr(a, row1, col1);
    printf("矩阵对角线的积\n");
    mul(a, row1, col1);
    printf("矩阵的1-范数\n");
    F1(a, row1, col1);
    printf("矩阵的无穷-范数\n");
    F2(a, row1, col1);
    printf("矩阵的Frobenius 范数\n");
    frobenius(a, row1, col1);
}

1.in

2 3

1 2 3

3 2 1

3 2

4 5

6 8

7 8

1.out

矩阵乘法

37 45

31 39

矩阵加法

Can't do this!

矩阵减法

Can't do this!

矩阵转置

1 3

2 2

3 1

矩阵的迹

Can't do this!

矩阵对角线的积

Can't do this!

矩阵的1-范数

6

矩阵的无穷-范数

4

矩阵的Frobenius 范数

5.2915

2.in

3 3

1 2 3

4 5 6

7 8 9

3 3

1 2 3

3 2 1

2 1 3

2.out

矩阵乘法

13 9 14

31 24 35

49 39 56

矩阵加法

2 4 6

7 7 7

9 9 12

矩阵减法

0 0 0

1 3 5

5 7 6

矩阵转置

1 4 7

2 5 8

3 6 9

矩阵的迹

15

矩阵对角线的积

45

矩阵的1-范数

24

矩阵的无穷-范数

18

矩阵的Frobenius 范数

16.8819

三、空间优化

对于任意的一个二维数组,我们可以采用一个一维数组进行定义,那么这个一维数组的大小即为row*col,我们可以在mian()函数中将其定义为。

double a[(row + 1) * (col + 1)];

此时,如果对于一个矩阵A我们要访问它的元素A_{i,j}可以通过数组a[i * (j - 1) + j] 来访问对应元素。
同时,在进行值传递的时候,可以略去空间声明,节省空间。下面给出其在矩阵乘法中的实现。

void dot(double a[], double b[], int row1, int col1, int row2, int col2)//矩阵乘法
{
    double c[(row1 + 1)*(col2 + 1)];
    memset(c, 0, sizeof(c));
    if (row2 != col1)
    {
        cout << "Can't do this!" << endl;
        cout << endl;
        return;
    }
    for (int i = 1; i <= row1; i++)
    {
        for (int j = 1; j <= col2; j++)
            for (int k = 1; k <= col1; k++)
                c[i * (j - 1) + j] += a[i * (k - 1) + k] * b[k * (j - 1) + j]; //答案存储矩阵
    }
    for (int i = 1; i <= row1; i++)//输出
    {
        for (int j = 1; j <= col2; j++)
            cout << c[i * (j - 1) + j] << ' ';
        cout << endl;
    }
    cout << endl;
}

对于其他部分的函数,转化方式类似,在此不进行赘述。


「雪霁融雾月,冰消凝夜雨」