矩阵与矩阵论

为啥要学矩阵？

       因为机器学习的大部分运算都是在矩阵上进行的，所以要对矩阵的各种运算比较熟，才能无障碍的理解和思考。本科线代的教学目标还是建立线性代数的各种基本概念，并不完全以矩阵的各种运算为主，比如矩阵求导，就完全没涉及。

      在学习信号处理、通信、模式识别的过程中，矩阵分析在科学研究中起到的重要作用。很多新的理论、方法和技术的诞生与发展就是矩阵理论和线性代数应用和推广的结果。如果学不好真的是寸步难行。其实真实点说就是，难以理解+难以研究+难以写论文。

     今天我们就来探讨一下矩阵的代码实现，并用两个小案例来分享一下。

2021数学建模A题 ——“FAST ”主动反射面的形状调节

先了解一下题目

      中国天眼——500 米口径球面射电望远镜（ Five hundred meter Aperture Spherical radioTelescope ，简称 FAST ），是我国具有自主知识产权 的目 前世界上单口径最大、灵敏 度最高 的射电望远镜。它的落成启用，对我国在科学沿实现重大原创突破、加快创新驱动发展具有重要意义。

  原题链接：https://www.aliyundrive.com/s/ekdxsvMc4bp

     今天主要分享一下这个第二题最后一问的代码思路（不讲数学的部分），俗话说的好啊，思维量越大，代码量越少。几何模型解法就是这样，代码在30行内就能搞定，但你的集合模型必须建立的可圈可点。

第二问：当待观测天体 位于 =36.795°,= 78.169°时， 确定理想抛物面。 建立反射面板调节模型，调节 相关促动器的伸缩量使反射面尽量贴近该理想抛物面。将理想抛物面的顶点坐标，以及调节后反射面 300 米口径内的主索节点编号、位置坐标 、各促动器的伸缩量等结果按照规定的格式（见附件 4 ）保存在 result.xlsx” 文件中。

     首先呢这是题物理题，最后的做法好像是有建立几何模型和算法方法，我的队伍呢没想到什么好的算法模型，所以用的是几何模型来解题。几何模型怎么和代码扯上关系呢？这就要回到我们今天的主题矩阵了。这时聪明的同学肯定已经想到了变换矩阵。

     概念：变换矩阵是数学线性代数中的一个概念。在线性代数中，线性变换能够用矩阵表示。如果T是一个把Rn映射到Rm的线性变换，且x是一个具有n个元素的列向量 ，那么我们把m×n的矩阵A，称为T的变换矩阵。

     最为常用的几何变换都是线性变换，这包括旋转、缩放、切变、反射以及正投影。在二维空间中，线性变换可以用 2×2 的变换矩阵表示。在本题的三维空间中，变换矩阵同样适用。

     简化可得解题思路：

1. 给原数据排序后，通过题目给出的"附件二.csv"（已知上下端点三维空间坐标）先求出促动器原长。
2. 将得出的几何模型简化成三维变换矩阵，使用Python读取csv数据文件进行变换运算，得到当前状态下理想抛物面对应的上端点新坐标，这即为题目要求中“主索节点调节后的位置坐标”。
3. 使用上端点的新三维坐标与促动器下端点（固定在地上，所以是不变的）的坐标计算出每个促动器新的长度。
4. 新旧长度相减即为题目要求“各促动器的伸缩量”。
5. 将运算得到的数据按照格式回填至excel表中提交，至此第二题结束。

     根据思路上代码：

这次用到的是numpy库，当然也有其他方法在py中表示矩阵。

先上三维坐标计算器计算原长：

import csv
def distance(one,two) :#原理就是三维空间的距离计算公式（根号下对应坐标差的平方）
    [x1,y1,z1] = one  # 第一个点的坐标
    [x2,y2,z2] = two  # 第二个点的坐标

    return (((x2-x1)**2)+((y2-y1)**2)+((z2-z1)**2))**(1/2)

with open('原长.csv', 'r') as f:
    reader = csv.reader(f)
    for row in reader:
        o = list(map(eval, row[2:5]))#数据切片根据题目给的格式设置
        t = list(map(eval, row[5:8]))
        print(distance(o,t))#输出每个单独促动器的原长，当然你也可以用xlwt库直接导入

运用变换矩阵得出新的坐标：

这里的变换矩阵不一定是最优解哈，四天三夜我队友给了我好多矩阵，随便挑了一个嘻嘻。

import csv
import numpy as np

temp = []
with open('附件2.csv', 'r') as f:
    reader = csv.reader(f)
    for row in reader:
        s1 = 0
        s2 = 0
        s3 = 0

        t = list(map(eval, row[5:8]))
        j = [t[0]]
        k = [t[1]]
        l = [t[2]]
        temp = [j,k,l]
        #print(temp)
        a = np.array([[0.1642,-0.97876, -0.1228], [0.7838,0.2050, - 0.5862], [0.5990, 0, 0.8008]])
#一个变换矩阵（随机挑的，不一定对哈），可以看出里面每个[]对应的是矩阵的每一行，数值是近似值
        c = np.array(temp)
        f1 = a * c #3*3的矩阵乘以3*1的矩阵，用的是叉乘
'''
这里提一嘴点乘是dot方法
d1=a@c
d2=tf.matmul(a,c)
d3=np.dot(a,c)
'''
#最核心的一行左乘变换矩阵得到新坐标，结果应该是一个3*1的矩阵，但实际测试却还是3*3的矩阵，是9个值。
#三个值一组的和为一个坐标，所以我加了一步求和，变为了三个正确的x,y,z的坐标。
        #print(f1)
        s1 = sum(f1[0])#s1代表坐标的第一列，后面以此类推
        s2 = sum(f1[1])
        s3 = sum(f1[2])
        print(s3)#这里把x,y,z的每一列坐标分开输出了，是为了方便我复制到excel，我懒得写导入的代码了QAQ

得出新的坐标之后代回第一段代码中，即可求出每一段新的长度，然后对应相减就可以得出各促动器的伸缩量了。

咱们第一个矩阵的代码案例也就分享结束了，接下来就是第二个，大名鼎鼎的兔子数题目。

【AI达人养成营】作业1 斐波那契数列——矩阵快速幂解法

     斐波那契数列（Fibonacci sequence），又称黄金分割数列，因数学家莱昂纳多·斐波那契（Leonardo Fibonacci）以兔子繁殖为例子而引入，故又称为“兔子数列”，指的是这样一个数列：0、1、1、2、3、5、8、13、21、34、……在数学上，斐波那契数列以如下被以递推的方法定义：F(0)=0，F(1)=1, F(n)=F(n - 1)+F(n - 2)（n ≥ 2，n ∈ N*）在现代物理、准晶体结构、化学等领域，斐波纳契数列都有直接的应用，为此，美国数学会从 1963 年起出版了以《斐波纳契数列季刊》为名的一份数学杂志，用于专门刊载这方面的研究成果。

先看题目：

好，既然是F1=F2=1，那么这就是个从1开始的数列，1，1，2，3，5，8……

只要注意第一个不是0就好了呗，那不是几行就结束了嘛。

n = int(input())
def fib(n):
    # 给递归一个出口  第一位和第二位都是1
    if n == 1 or n == 2:
        return 1
    else:
        # 从第三位开始  返回上一个数加上上一个数
        return fib(n-1) + fib(n-2)
print(fib(n)%10007)

定睛一看呢，不对啊，这题坏的很，诶，想要拿满分得满足时间限制，而递归肯定会超时。

重新读题：

那为了不超时，只能请出我们的主角，矩阵快速幂了。

假设我要求A的9次方：

A^9

= A*A*A*A*A*A*A*A*A  【一个一个乘，要乘9次】

= A*(A*A)*(A*A)*(A*A)*(A*A) = A*(A^2)^4 【A平方后，再四次方，还要乘上剩下的一个A，要乘6次】

= A*((A^2)^2)^2【A平方后，再平方，再平方，还要乘上剩下的一个A，要乘4次】

1000000次幂，暴力要乘1000000次，快速幂就只要（log2底1000000的对数） 次，约等于20次，这大大降低了程序的运行时间。

那么对于这道题呢：

核心思想就这一个矩阵：

它的时间复杂度为 O(log n)，所以可以完成任务！

具体实现代码如下（无需numpy，用循环递推）：

n = int(input())
def mul(a, b):# 首先定义二阶矩阵乘法运算
    c = [[0, 0], [0, 0]]# 定义一个空的二阶矩阵，存储结果
    for i in range(2):# row
        for j in range(2):# col
            for k in range(2):# 新二阶矩阵的值计算
                c[i][j] += a[i][k] * b[k][j]
    return c

def fib(n):
    if n <= 1: return max(n, 0)
    temp = [[1, 0], [0, 1]]# 单位矩阵，等价于1
    A = [[1, 1], [1, 0]]# 得出的base矩阵，带入进行运算
    while n:
        if n & 1: temp = mul(temp, A)# 如果n是奇数，或者直到n=1停止
        A = mul(A, A)#对A进行快速幂
        n >>= 1 #整除2，向下取整，这里使用数据右移更快
    return temp[0][1]

result = fib(n)%10007
print(result)

这题还有公式法、暴力求解、递归、生成器等方法大家可以自行探索。