-
什么是矩阵
2019-03-19 16:55:24矩阵,在数学上,矩阵是指纵横排列的二维数据表格,最早来自于方程组的系数及常数所构成的方阵。矩阵是高等代数学中的常见工具,也常见于统计分析等应用数学学科中。 长这个样子: 失量也可以转为矩阵,可以看成nX1...什么是矩阵
矩阵,在数学上,矩阵是指纵横排列的二维数据表格,最早来自于方程组的系数及常数所构成的方阵。矩阵是高等代数学中的常见工具,也常见于统计分析等应用数学学科中。
长这个样子:
矢量也可以转为矩阵,可以看成nX1的行矩阵,或1Xn的矩阵。
矩阵列的运行比较复杂,下面就来一一探讨。矩阵和标量的乘法
直接标量与各个分量相乘即可,不多废话了…同时kM=Mk即,谁在哪边都一样。
矩阵与矩阵的乘法
它会得到一个新的矩阵,而且维度与这两个矩阵有关系。
如A为4X3矩阵,B为3X6矩阵那么 AB维度就是4X6。
左矩阵的列数必须与右矩阵的行数想同,否则不能相乘。
矩阵不满足交换律:AB!=BA
满足结合律:(AB)C=A(BC) 甚至可以扩展至 ABCDE=((A(BC))D)E=(AB)(CD)E方阵
方块矩阵,即行列数相同的矩阵。有一些运算和性质是只有方阵有具有,如对角元素
对角矩阵:
单位矩阵(I):
单位矩阵乘完还等于原本的矩阵,设I为转:
MI=IM=M
转置矩阵(Mt)
对原矩阵的一种运算,即行变列,列变行。可以记作Mt
性制一:转两次就转回来了:
(Mt)t=M
性制二:矩阵串接转置,等于反射串接各矩阵
(AB)t=BtAt逆矩阵(M-1)
这应该是这里最复杂的一种操作了。不是所有矩阵都有逆矩阵,它必须是一个方阵。
给定M-1来表示。最重要的特性就是M和M-1相乘会得到一个单位矩阵。也就是说:
MM-1=M-1M=I并非所有有对应的逆矩阵,如果一个矩阵有对应的逆矩阵则这个矩阵称为是可逆的,否则称为不可逆的。
如果一个矩阵行列式不为0,那么它就是可逆的。性质一:逆矩阵的逆矩阵就是它本身
(M-1)-1=M性质二:单位矩阵的逆矩阵就是它本身
I-1=I性制三:转置矩阵的逆矩阵是逆矩阵的转置
(Mt)-1=(M1)t性质四:矩阵串接相乘后的逆矩阵等于反向串接各个矩阵的逆矩阵
(ABCD)-1=D-1C-1B-1A-1性质五:允许我们还原这个变换
M-1(Mv)=(M-1M)v=Iv=v正交矩阵
方正M和它的转置矩阵乘积为单位矩阵的话,它就是一个正交矩阵,即:
MMt=MtM=I正交矩阵的逆矩阵和转置矩阵是一样的
Mt=M-1三维变换中我们经常会需要作用逆矩阵来求解反射的变换。而逆矩阵的求解往往计算量很大,但转置矩阵就非常容易。
矩阵的每一行,即c1、c2、c3的是单位矢量,由于其相互垂直只有与自己点乘才能得到1,其他为0.矩阵与矢量相乘
我们需要把矢量先转成行矩阵或是列矩阵,但要满足矩阵相乘的条件。通常我们使用右乘。
-
矩阵乘法的常数优化
2015-04-16 10:03:44矩阵乘法的常数优化 philipsweng 虽然说作为键盘科学家,我们更应该关心程序的时间复杂度。但是一个写的不好的程序可能在实际运行会跟时间复杂度更差的程序差不了多少。...这是为什么??? 注意到第二个程序的k矩阵乘法的常数优化 philipsweng
虽然说作为键盘科学家,我们更应该关心程序的时间复杂度。但是一个写的不好的程序可能在实际运行会跟时间复杂度更差的程序差不了多少。我们我们也应该注意程序的常数优化。
对于矩阵乘法来讲。我们实际上可以比较这两种打法。实际上第二种打法在绝大多数情况下效率约为第一种的两倍。
这是为什么???
注意到第二个程序的k放在了第二层。
C[i][j] = a[i][k] * b[k][j]
那么我们枚举的顺序就使得a[i]数据顺序枚举以及b[k]顺序枚举。
那么我们数组的寻址速度会增加许多。
然后对于取模操作,我们并不需要每次都取模,因为取模太缓慢。我们可以界定一个范围再取模(慎用。容易溢出) -
随机变量和常数的协方差等于_神奇又好玩的协方差矩阵
2021-01-08 15:32:00今天我们不讲雷达,我们...所以今天我们就来把基础打好,学一学这个神奇又好玩的数学工具——协方差矩阵一维数据的分析方法——方差方差是在概率论和统计方差衡量随机变量或一组数据时离散程度的度量,样本方差的无...今天我们不讲雷达,我们来讲讲协方差矩阵。为什么呢?因为如果我直接的去说超分辨估计的东西,放在一篇文章里面太长了,我自己写着都累,别人看着也累。同时协方差矩阵这个东西又比较基础,不仅雷达要用,其他各个工程领域也喜欢这个数学工具。所以今天我们就来把基础打好,学一学这个神奇又好玩的数学工具——协方差矩阵
一维数据的分析方法——方差
方差是在概率论和统计方差衡量随机变量或一组数据时离散程度的度量,样本方差的无偏估计公式为:
sigma = 1; x = sigma .* randn(100000,1); figure(1) ksdensity(x);grid on;title('概率分布 sigma = 1')
方差能很好的解释了一维数据的分布特性,但是当我们开始分析二维数据的时候,完球了!
新问题——二维到N维数据
我们来构造一组二维数据:
x = 4.*randn(1,1000); y = 1.*randn(1,1000); D = [x;y]; plot(D(1,:),D(2,:),'.'); grid on; axis equal;
我们构造的这组数据,一个是沿着X轴标准差为4和以及沿着Y轴标准差为1的数据,画图看到是这样的,但是这样子其实过于简单,大多数真实情况是比较复杂的,这种根本解决不了问题,为了增加点难度,我们引入了多元正态分布与线性变换的概念。
更真实的模拟——多元正态分布与线性变换
我们再搞点事情,用旋转矩阵把数据给他搅一下,实际上是我们在对原始数据做了一个变换
,矩阵
变换矩阵(transformation matrix),为了将图中的点经过线性变换得到我们想要的,我们需要构造这个矩阵被称为
例如有尺度矩阵
或者旋转矩阵
其中
是旋转的度数,那么在这个例子中我们用旋转矩阵来折腾数据,随手写个旋转矩阵函数
function out =rotate_2D(thetad) out =[cosd(thetad)-sind(thetad);sind(thetad)cosd(thetad)]; end
对原始数据进行旋转
N = 1000; x = 4.*randn(1,N); y = 1.*randn(1,N); D = [x;y]; new_D = rotate_2d(40)*D; plot(new_D(1,:),new_D(2,:),'.');grid on;axis equal
或者可以这么写,用矩阵的处理方法
N = 1000; x = 1.*randn(1,N); y = 1.*randn(1,N); D = [x;y]; T = rotate_2d(40)*[4 0;0 1]; new_D = T*D; plot(new_D(1,:),new_D(2,:),'.');grid on;axis equal
上面两种写法是等价的,都可以,不过为了理解方便,我用尺度矩阵和旋转矩阵继续往下整(挖坑了)。那么这个4和1就是标准差的缩放了,因此变换矩阵可以表示为
如果式中
,我们会看到下图:
这下难度够了。二维数据,而且并没有像之前那样子是按照X和Y轴分布的,比较接近真实情况,构造数据的问题解决了,那么下来怎么分析呢?
新工具——协方差/协方差矩阵
看到这么一种数据,如果按照我们如果按照一维数据的分析思路方法,计算在X方向上的方差
和Y方向上的方差
两个轴分解产生了正相关。怎么办?,肯定是会出问题。因为图形明显的告诉你,单纯的从水平方向和垂直方向去分析并不能对这个斜着分布的关系做出合理的解释,并且x值的增加同时y值也在增加,沿着这
这时候我们就要扩展方差的概念了,为了捕捉这个特征,一帮闲的没事的人们就开始整活——协方差出现了:
对于二维的数据,我们可以得到的有这四个东西
我们用矩阵来表示,这个矩阵称为协方差矩阵:
由于X与Y的相关性和Y与X的相关性是一样的,因此
协方差矩阵始终是一个对称阵,其对角线上的是方差,非对角线上的是协方差。二维正态分布的数据由它的均值和 2x2 的协方差矩阵就可以表示,我们将刚才的数据协方差矩阵打印出来:
Rxx = new_D *new_D './(N-1)
随着数据形态的变化,协方差矩阵的值也跟着在变,大家可以在代码里面修改一些参数试试,我就不上图了,那么问题又来了,协方差矩阵有了,我们要怎么分析,才能较好的还原数据的原始特性呢?
协方差矩阵的特征分解
从前面的部分可以知道,协方差矩阵定义了数据的离散程度(方差)和空间分布走向(协方差)。因此,我们可以认为,指向最大离散度方向的向量和其值(该值等于这个方向上的方差),表示协方差矩阵。
假设前面的数据是
,实际上我们刚才的操作其实是对数据
进行线性变换,可以得到:
其中
是含有旋转矩阵
和缩放矩阵
的变换矩阵:
矩阵
和
的定义我们再贴出来复习一下:
接下来我们就要讨论协方差矩阵
和线性变换
之间的关系。
我们来回忆一下刚才的数据构造流程,先从一组无缩放的数据开始,也就是标准差为1的无旋转,服从正态分布的数据,我们可以称之为"白数据”。
白数据的协方差矩阵,等于单位阵,因此方差和标准差为1,协方差为0:
现在,对
数据的标准差缩放到4:
变换后的数据
如图所示:
的协方差矩阵
为:
故
的协方差矩阵
,与作用在原始数据上的线性变换矩阵
的关系为
,并且:
但是这个式子是否适用于旋转变换呢?我们将协方差矩阵分解为旋转和缩放矩阵的乘积。
协方差矩阵可以用特征向量和特征值表示,在二维情况下,有2个特征向量和2个特征值。通过矩阵来定义两个等式的关系:
式中
是由
的各特征向量构成的矩阵,
为相应特征向量对应的非0特征值组成的对角矩阵,上式可以写为我们常见的形式:
这个式子我们称为协方差矩阵的特征分解。由于特征向量表示了最大方差的方向,特征值表示了该方向上方差的值,也就是说,
实际上就是旋转矩阵,而
协方差矩阵可以进一步表示为:就是缩放矩阵,
式中
为旋转矩阵,
为缩放矩阵。因为
为对角缩放矩阵,所以
,并且由于
为正交矩阵,所以
。那么协方差矩阵最后的形式可以写为:
也就是说,协方差矩阵特征分解后的特征向量,总是指向数据最大方差的方向,并规定了其走向。而且由于旋转矩阵的正交性,次最大特征值对应的特征向量,总是与最大特征值的特征向量正交。
我们把代码补完,验证一下:
N = 1000; x = 1.*randn(1,N); y = 1.*randn(1,N); D = [x;y]; T = rotate_2d(40)*[4 0;0 1]; new_D = T*D; Rxx = new_D*new_D'./(N-1); [V,D] = eig(Rxx) figure(1) plot(new_D(1,:),new_D(2,:),'.');grid on;axis equal;hold on; quiver(0,0,V(1,1).*sqrt(D(1,1)),V(2,1).*sqrt(D(1,1)),'LineWidth',3) quiver(0,0,V(1,2)*sqrt(D(2,2)),V(2,2)*sqrt(D(2,2)),'LineWidth',3) legend('数据分布','特征向量1','特征向量2') hold off
我们看到特征向量其实就是原来的旋转矩阵,而特征值的结果也符合原来数据的方差,即1和16,即为原来的标准差1和4。
那么更高维度的数据是否也符合要求呢?
我们构造一个三维数据试一试:
旋转函数:
function yaw = rotate_yd(theta) yaw = [cosd(theta) 0 sind(theta);0 1 0;-sind(theta) 0 cosd(theta)]; end
代码部分:
clear;clc;format longG %-------------------------------------------------------------------------- % 一组标准的概率分布,构造一个二维数据 %-------------------------------------------------------------------------- N = 2000; x = 5.*randn(1,N); y = 1.*randn(1,N); z = 3.*randn(1,N); A = [x;y;z]; A = rotate_zd(0)*rotate_yd(20)*rotate_xd(0)*A; Rxx = A*A'./(N-1); [V,D] = eig(Rxx); figure(1) plot3(A(1,:),A(2,:),A(3,:),'.');grid on;axis equal;hold on; quiver3(0,0,0,V(1,1).*sqrt(D(1,1)).*5,V(2,1).*sqrt(D(1,1)).*5,V(3,1).*sqrt(D(1,1)).*5,'LineWidth',2) quiver3(0,0,0,V(1,2)*sqrt(D(2,2)).*5,V(2,2)*sqrt(D(2,2)).*5,V(3,2)*sqrt(D(2,2)).*5,'LineWidth',2) quiver3(0,0,0,V(1,3)*sqrt(D(3,3)).*5,V(2,3)*sqrt(D(3,3)).*5,V(3,3)*sqrt(D(3,3)).*5,'LineWidth',2) hold off
适当的缩放了一下特征向量的长度,我乘了个5进去,为了画图好看,我们看看图形:
总结
数据的协方差矩阵,直接与数据的线性变换有关,该线性变换完全由特征向量和特征值定义,特征向量表示了矩阵的旋转,特征值对应各个维度上的缩放因子的平方。
参考:
Xinyu Chen:如何直观地理解「协方差矩阵」?zhuanlan.zhihu.comA geometric interpretation of the covariance matrixwww.visiondummy.com翻译:协方差矩阵的几何解释njuferret.github.io作为2020年的最后一天,以这篇文章完满结束收尾。我是从今年才开始决定写文章,想把自己的一些学习心得更新上来,但是自己又比较拖延,速度比较慢,希望大家能够见谅,写一篇文章下来需要好多天,每篇都保证较好的图,文字和公式说明,所以确实是一件花心思的事情。
那么希望在2021年里,能够继续保持下去,持续的把自己在工作中学习中的一些理解分享给大家,能够帮助到你们就是最大的快乐。
最后,有不足之处欢迎指正,评论区不能及时回复,望见谅,转载请注明知乎来源,谢谢~
-
矩阵乘法
2017-09-21 15:01:43矩阵是什么 矩阵加减法 矩阵乘法 矩阵乘法的运算定律 模板 矩阵是什么? 就是矩阵 矩阵加减法? 矩阵加矩阵,就将对应的加上。 矩阵加常数,就将每一个元素都加上这个常数。 减法同理。 矩阵乘法 ...矩阵乘法是个很高端的东西。
注意事项: 以下的讲述下标都从0开始矩阵是什么?
不解释
矩阵加减法?
矩阵加矩阵,就将对应的加上。
矩阵加常数,就将每一个元素都加上这个常数。
减法同理。矩阵乘法
一张好图,在这里发现的,方便理解矩阵乘法。
一个m∗n 的的A矩阵,和一个n∗p 的B矩阵相乘,将得到一个m∗p 的矩阵C
C(i,j)=∑k=1nA(i,k)∗B(k,j)
可以简单地理解为,A中i行的元素,与B中j列的元素,对应相乘得到的和。矩阵乘法的运算定律
(1) 不满足交换律
(2) 满足结合律:(AB)C=A(BC)
(3) 满足分配律:(A+B)C=AC+BC A(B+C)=AB+AC模板
template <typename T,int M,int N> struct Matrix { T mat[M][N]; inline Matrix(){memset(mat,0,sizeof mat);} inline T* operator[](int i){return mat[i];} }; template <typename T,int M,int N,int P> inline Matrix<T,M,P> operator*(const Matrix<T,M,N>& x,const Matrix<T,N,P>& y) { Matrix<T,M,P> ret; int i,j,k; for (i=0;i<M;++i) for (j=0;j<P;++j) for (k=0;k<N;++k) ret.mat[i][j]+=x.mat[i][k]*y.mat[k][j]; return ret; } template <typename T,typename T_> inline T pow(T x,T_ y) { T ret=x; --y;//这两句话可以忽略初值问题,但y<=0时就悲剧了 while (y) { if (y&1) ret=ret*x; y>>=1; x=x*x; } return ret; }
矩阵乘法的应用
矩阵乘法可以优化时间。
将某些操作转化为矩阵乘法的形式,就可以用快速幂减少时间。因为矩阵乘法满足结合律。
例题一 斐波那契数列
描述
题目背景
大家都知道,斐波那契数列是满足如下性质的一个数列:
• f(1) = 1
• f(2) = 1
• f(n) = f(n-1) + f(n-2) (n ≥ 2 且 n 为整数)题目描述
请你求出 f(n) mod 1000000007 的值。
输入格式:
·第 1 行:一个整数 n
输出格式:
第 1 行: f(n) mod 1000000007 的值
输入样例#1:
5
输出样例#1:
5
输入样例#2:
10
输出样例#2:
55
说明
对于 60% 的数据: n ≤ 92
对于 100% 的数据: n在long long(INT64)范围内。解法
这是一道裸斐波拉契数列。传统递推是O(N)的,显然会炸。
斐波拉契数列有个通项,但我们在这里用矩阵乘法解决。
我们知道f(n)是由f(n-1)和f(n-2)推过来的,不妨设一个1*2的矩阵
A=[f(n−2)f(n−1)]
现在我们要通过A推出B
B=[f(n−1)f(n)]=[f(n−1)f(n−2)+f(n−1)]
我们要求出一个2*2的 转移矩阵 T,满足
A∗T=B
即
[f(n−2)f(n−1)]∗T=[f(n−1)f(n−2)+f(n−1)]
我们知道,B[0][0]=A[0][1] B[0][1]=A[0][0]+A[0][1]
对于T[0][0],A[0][0]对B[0][0]没有贡献,所以T[0][0]=0
对于T[1][0],A[0][1]对B[0][0]有贡献,B[0][0]=A[0][1]*1,所以T[1][0]=1
对于T[0][1],A[0][0]对B[0][1]有贡献,B[0][1]=A[0][0]*1+A[0][1]*1,所以T[0][1]=1
对于T[1][1],A[1][1]对B[0][1]有贡献,B[1][1]=A[0][0]*1+A[0][1]*1,所以T[1][1]=1
综上所述,
T=[0111]
显然这个式子是正确的。
[f(n−2)f(n−1)]∗[0111]=[f(n−1)f(n−2)+f(n−1)]
求出这个矩阵后,就可以愉快地快速幂了。
时间复杂度O(lgN)代码
#include <cstdio> #include <cstring> #include <algorithm> using namespace std; #define mod 1000000007 template <typename T,int M,int N> struct Matrix { T mat[M][N]; inline Matrix(){memset(mat,0,sizeof mat);} inline T* operator[](int i){return mat[i];} }; template <typename T,int M,int N,int P> inline Matrix<T,M,P> operator*(const Matrix<T,M,N>& x,const Matrix<T,N,P>& y) { Matrix<T,M,P> ret; int i,j,k; for (i=0;i<M;++i) for (j=0;j<P;++j) for (k=0;k<N;++k) (ret.mat[i][j]+=x.mat[i][k]*y.mat[k][j])%=mod; return ret; } template <typename T,typename T_> inline T pow(T x,T_ y) { T ret=x; --y; while (y) { if (y&1) ret=ret*x; y>>=1; x=x*x; } return ret; } int main() { long long n; scanf("%lld",&n); if (n==1) putchar('1'); else { Matrix<long long,2,2> tmp; tmp[0][1]=tmp[1][0]=tmp[1][1]=1; tmp=pow(tmp,n-1); Matrix<long long,1,2> a; a[0][1]=1;//a=[f(0) f(1)] a=a*tmp;//[f(0) f(1)]*(T^(n-1))=[f(n-1) f(n)] printf("%lld\n",a[0][1]); } }
例题二 【NOIP2013模拟联考14】图形变换
Description
翔翔最近接到一个任务,要把一个图形做大量的变换操作,翔翔实在是操作得手软,决定写个程序来执行变换操作。
翔翔目前接到的任务是,对一个由n个点组成的图形连续作平移、缩放、旋转变换。相关操作定义如下:
Trans(dx,dy) 表示平移图形,即把图形上所有的点的横纵坐标分别加上dx和dy;
Scale(sx,sy) 表示缩放图形,即把图形上所有点的横纵坐标分别乘以sx和sy;
Rotate(θ,x0,y0) 表示旋转图形,即把图形上所有点的坐标绕(x0,y0)顺时针旋转θ角度
由于某些操作会重复运行多次,翔翔还定义了循环指令:
Loop(m)
…
End
表示把Loop和对应End之间的操作循环执行m次,循环可以嵌套。Input
第一行一个整数n(n<=100)表示图形由n个点组成;
接下来n行,每行空格隔开两个实数xi,yi表示点的坐标;
接下来一直到文件结束,每行一条操作指令。保证指令格式合法,无多余空格。Output
输出有n行,每行两个空格隔开实数xi,yi表示对应输入的点变换后的坐标。
本题采用Special Judge判断,只要你输出的数值与标准答案误差不能超过1即可。Sample Input
3
0.5 0
2.5 2
-4.5 1
Trans(1.5,-1)
Loop(2)
Trans(1,1)
Loop(2)
Rotate(90,0,0)
End
Scale(2,3)
EndSample Output
10.0000 -3.0000
18.0000 15.0000
-10.0000 6.0000Data Constraint
保证操作中坐标值不会超过double范围,输出不会超过int范围;
指令总共不超过1000行;
对于所有的数据,所有循环指令中m<=1000000;
对于60%的数据,所有循环指令中m<=1000;
对于30%的数据不含嵌套循环。Hint
【友情提醒】
pi的值最好用系统的值。C++的定义为:#define Pi M_PI
Pascal就是直接为:pi
不要自己定义避免因为pi带来的误差。解法
旋转公式:
绕(a,b)逆时针转θx′=(x−a)∗cosθ−(y−b)∗sinθ+ay′=(x−a)∗sinθ+(y−b)∗cosθ+b
暴力铁定超时。
先考虑使用2*2的转移矩阵T,但是只能对付乘法,加法和旋转就要用矩阵加法了。矩阵乘法和矩阵加法混在一起很麻烦(矩阵套矩阵?)。
我们可以新增一个常数项
A=[xy1]
现在要求出T1、T2、T3满足
[xy1]∗T1=[x+ay+b1][xy1]∗T2=[x∗ay∗b1][xy1]∗T3=[(x−a)∗cosθ−(y−b)∗sinθ+a(x−a)∗sinθ+(y−b)∗cosθ+b1]
这样加、乘就特别容易,将旋转公式用乘法分配律拆开,也可以得出转移矩阵
T1=⎡⎣⎢10a01b001⎤⎦⎥T2=⎡⎣⎢a000b0001⎤⎦⎥T3=⎡⎣⎢cosθ−sinθa−a∗cosθ+b∗sinθsinθcosθb−a∗sinθ−b∗cosθ001⎤⎦⎥
可以代进去看看,是成立的。
于是我们可以将所有操作乘起来(有循环时递归,算出这个循环里一次的乘积后快速幂),最后在将每个坐标乘这个积,就是答案。代码
英语不好,将theta打成xida,懒得改,注意一下就好了。
#include <cstdio> #include <cstring> #include <cmath> #include <cassert> #include <algorithm> using namespace std; #define I_O(x) freopen(""#x".in","r",stdin);freopen(""#x".out","w",stdout) #define Pi 3.14159265358979323846 template <typename T,int M,int N> struct Matrix { T mat[M][N]; inline Matrix(){memset(mat,0,sizeof mat);} inline T* operator[](int i){return mat[i];} }; template <typename T,int M,int N,int P> inline Matrix<T,M,P> operator*(const Matrix<T,M,N>& x,const Matrix<T,N,P>& y) { Matrix<T,M,P> ret; int i,j,k; for (i=0;i<M;++i) for (j=0;j<P;++j) for (k=0;k<N;++k) ret.mat[i][j]+=x.mat[i][k]*y.mat[k][j]; return ret; } template <typename T,typename T_> inline T pow(T x,T_ y) { T ret=x; --y; while (y) { if (y&1) ret=ret*x; y>>=1; x=x*x; } return ret; } int n; Matrix<double,1,3> d[100]; char cz[101]; Matrix<double,3,3> dfs(int); int main() { I_O(transform); scanf("%d",&n); int i; double x,y; for (i=0;i<n;++i) { scanf("%lf%lf",&d[i][0][0],&d[i][0][1]); d[i][0][2]=1; } Matrix<double,3,3> a,tmp1,tmp2,tmp3;//分三个变量是为避免多次赋值 a[0][0]=a[1][1]=a[2][2]=1;//a的初值 tmp1[0][0]=tmp1[1][1]=tmp1[2][2]=1; tmp2[2][2]=1; tmp3[2][2]=1; double xida,tmp_sin,tmp_cos; int t; while (scanf("%s",cz)==1) { switch (*cz) { case 'T': // 1 0 0 // 0 1 0 // a b 1 sscanf(cz,"Trans(%lf,%lf",&tmp1[2][0],&tmp1[2][1]); a=a*tmp1; break; case 'S': // a 0 0 // 0 b 0 // 0 0 1 sscanf(cz,"Scale(%lf,%lf",&tmp2[0][0],&tmp2[1][1]); a=a*tmp2; break; case 'R': sscanf(cz,"Rotate(%lf,%lf,%lf",&xida,&x,&y); xida=-xida/180*Pi;//顺时针角度转成逆时针弧度 // cos sin 0 // -sin cos 0 // a-a*cos+b*sin b-a*sin-b*cos 1 tmp_sin=sin(xida); tmp_cos=cos(xida); tmp3[0][0]=tmp_cos; tmp3[0][1]=tmp_sin; tmp3[0][2]=0; tmp3[1][0]=-tmp_sin; tmp3[1][1]=tmp_cos; tmp3[1][2]=0; tmp3[2][0]=x-x*tmp_cos+y*tmp_sin; tmp3[2][1]=y-x*tmp_sin-y*tmp_cos; a=a*tmp3; break; default: sscanf(cz,"Loop(%d",&t); a=a*dfs(t); } } Matrix<double,1,3> tmp; for (i=0;i<n;++i) { tmp=d[i]*a; printf("%.4lf %.4lf\n",tmp[0][0],tmp[0][1]); } } Matrix<double,3,3> dfs(int t) { double x,y; Matrix<double,3,3> a,tmp1,tmp2,tmp3; a[0][0]=a[1][1]=a[2][2]=1; tmp1[0][0]=tmp1[1][1]=tmp1[2][2]=1; tmp2[2][2]=1; tmp3[2][2]=1; int times; double xida,tmp_sin,tmp_cos; for (scanf("%s",cz);*cz!='E';scanf("%s",cz)) { switch (*cz) { case 'T': sscanf(cz,"Trans(%lf,%lf",&tmp1[2][0],&tmp1[2][1]); // 1 0 0 // 0 1 0 // a b 1 a=a*tmp1; break; case 'S': sscanf(cz,"Scale(%lf,%lf",&tmp2[0][0],&tmp2[1][1]); // a 0 0 // 0 b 0 // 0 0 1 a=a*tmp2; break; case 'R': sscanf(cz,"Rotate(%lf,%lf,%lf",&xida,&x,&y); xida=-xida/180*Pi; // cos sin 0 // -sin cos 0 // a-a*cos+b*sin b-a*sin-b*cos 1 tmp_sin=sin(xida); tmp_cos=cos(xida); tmp3[0][0]=tmp_cos; tmp3[0][1]=tmp_sin; tmp3[0][2]=0; tmp3[1][0]=-tmp_sin; tmp3[1][1]=tmp_cos; tmp3[1][2]=0; tmp3[2][0]=x-x*tmp_cos+y*tmp_sin; tmp3[2][1]=y-x*tmp_sin-y*tmp_cos; a=a*tmp3; break; default: sscanf(cz,"Loop(%d",×); a=a*dfs(times); } } return pow(a,t); }
例题三 【SDOI2009】HH去散步
描述
题目描述
HH有个一成不变的习惯,喜欢饭后百步走。所谓百步走,就是散步,就是在一定的时间 内,走过一定的距离。 但是同时HH又是个喜欢变化的人,所以他不会立刻沿着刚刚走来的路走回。 又因为HH是个喜欢变化的人,所以他每天走过的路径都不完全一样,他想知道他究竟有多 少种散步的方法。
现在给你学校的地图(假设每条路的长度都是一样的都是1),问长度为t,从给定地 点A走到给定地点B共有多少条符合条件的路径输入格式:
第一行:五个整数N,M,t,A,B。其中N表示学校里的路口的个数,M表示学校里的 路的条数,t表示HH想要散步的距离,A表示散步的出发点,而B则表示散步的终点。
接下来M行,每行一组Ai,Bi,表示从路口Ai到路口Bi有一条路。数据保证Ai != Bi,但 不保证任意两个路口之间至多只有一条路相连接。 路口编号从0到N − 1。 同一行内所有数据均由一个空格隔开,行首行尾没有多余空格。没有多余空行。 答案模45989。输出格式:
一行,表示答案。
输入样例#1:
4 5 3 0 0
0 1
0 2
0 3
2 1
3 2输出样例#1:
4
说明
对于30%的数据,N ≤ 4,M ≤ 10,t ≤ 10。
对于100%的数据,N ≤ 50,M ≤ 60,t ≤ 2^30,0 ≤ A,B解法
先说一下大意(我等好久才弄懂)。
给你一张图,要从A走到B,不能走上一次走的路,要走t步,问方案数。
不能走上一次走的路,一次。
举个例子:
路径可以这样:0->1->2->3->1->0
所以我们可以想到DP
设f[i][j]表示第i步走了j这条边的方案数,这样在转移时,就可以方便地判断,避免走上一次走的路。
转移:f[i][bc]=∑ab≠bcf[i−1][ab]
但是,我们发现,每次的决策点都是一样的,很明显能直接用矩阵乘法来转移,快速幂即可。代码
#include <cstdio> #include <cstring> #include <cstdlib> #include <algorithm> using namespace std; template <typename T,int M,int N> struct Matrix { T mat[M][N]; inline Matrix(){memset(mat,0,sizeof mat);} inline T* operator[](int i){return mat[i];} }; template <typename T,int M,int N,int P> inline Matrix<T,M,P> operator*(const Matrix<T,M,N>& x,const Matrix<T,N,P>& y) { Matrix<T,M,P> ret; int i,j,k; for (i=0;i<M;++i) for (j=0;j<P;++j) for (k=0;k<N;++k) ret.mat[i][j]+=x.mat[i][k]*y.mat[k][j]; for (i=0;i<M;++i) for (j=0;j<P;++j) ret.mat[i][j]%=45989;//减少mod运算(mod运算很慢的) return ret; } template <typename T,typename T_> inline T pow(T x,T_ y) { T ret=x; --y; while (y) { if (y&1) ret=ret*x; y>>=1; x=x*x; } return ret; } int n,m,t,u,v; struct EDGE { int from,to; EDGE* las; EDGE* rev; } e[120]; EDGE* last[20]; Matrix<long long,1,120> f; Matrix<long long,120,120> T; int li[60];//用于记录连接终点的边 int nl; int main() { scanf("%d%d%d%d%d",&n,&m,&t,&u,&v); int i,j=-1,x,y; for (i=1;i<=m;++i) { scanf("%d%d",&x,&y); ++j; e[j]={x,y,last[x],e+j+1}; last[x]=e+j; if (y==v) li[nl++]=j; ++j; e[j]={y,x,last[y],e+j-1}; last[y]=e+j; if (x==v) li[nl++]=j; } EDGE* ei; for (ei=last[u];ei;ei=ei->las) f[0][int(ei-e)]=1; m<<=1; EDGE *ej,*end=e+m; for (ei=e;ei<end;++ei) for (ej=last[ei->to];ej;ej=ej->las) if (ei->rev!=ej) T[int(ei-e)][int(ej-e)]=1;//生成转移矩阵 f=f*pow(T,t-1); int ans=0; for (i=0;i<nl;++i) ans+=f[0][li[i]]; printf("%d\n",ans%45989); }
总结
- 遇到矩阵乘法的题,先将转移后的式子拆开,再帮它配转移矩阵。
- 如果用普通的方式配矩阵必须出现加法,那么可以尝试加常数项上去。
- 若有某些DP(还是递推?)的题目的决策点固定,那么可以在程序中生成转移矩阵,直接快速幂。
补充
矩阵乘法模板2.0
template <typename T,int M,int N> struct Matrix { int m,n; T mat[M][N]; inline Matrix(){m=M;n=N;memset(mat,0,sizeof mat);} inline void SetUpSize(int _m,int _n){m=_m;n=_n;} inline T* operator[](int i){return mat[i];} }; template <typename T,int M,int N,int P> inline Matrix<T,M,P> operator*(const Matrix<T,M,N>& x,const Matrix<T,N,P>& y) { Matrix<T,M,P> ret; int i,j,k; for (i=0;i<x.m;++i) for (j=0;j<y.n;++j) for (k=0;k<x.n;++k) ret.mat[i][j]+=x.mat[i][k]*y.mat[k][j]; return ret; } template <typename T,typename T_> inline T pow(T x,T_ y) { T ret=x; --y;//这两句话可以忽略初值问题,但y<=0时就悲剧了 while (y) { if (y&1) ret=ret*x; y>>=1; x=x*x; } return ret; }
-
矩阵分析 第二章 lambda矩阵和Jordan标准型
2018-11-17 21:05:13什么是矩阵 矩阵元素是的多项式就是矩阵 矩阵的秩 行列式可以是的多项式但是能是零,就可以求 矩阵的逆 首先,求它的行列式,只有行列式为常数且非0;这个矩阵才有逆。 然后求伴随矩阵,再除于这个行列式 ... -
吴恩达的机器学习--矩阵运算
2019-06-05 08:44:35在数学中,矩阵(Matrix)是一个按照长方阵列排列的复数或实数集合,最早来自于方程组的系数及常数所构成的方阵。 是高等代数学中的常见工具,也常见于统计分析等应用数学学科中。 矩阵:由数字组成的矩形队列 ... -
matdem矩阵维度必须一致_吴恩达的机器学习--矩阵运算
2020-12-31 08:40:38什么是矩阵?在数学中,矩阵(Matrix)是一个按照长方阵列排列的复数或实数集合,最早来自于方程组的系数及常数所构成的方阵。是高等代数学中的常见工具,也常见于统计分析等应用数学学科中。矩阵:由数字组成的矩形... -
opengl 旋转矩阵和纹理坐标相乘_游戏引擎养成《六》 数学库-矩阵
2020-12-22 16:50:14# 什么是矩阵 在数学中,矩阵(Matrix)是一个按照长方阵列排列的复数或实数集合 ,最早来自于方程组的系数及常数所构成的方阵。 这一概念由19世纪英国数学家凯利首先提出。矩阵是高等代数学中的常见工具,也常见于... -
【AR】计算摄像头相对于探测到的标识的转换矩阵中的矩阵
2018-05-19 09:50:34在ARToolkit矩阵中的第四部就是计算摄像头相对于探测到的标识的转换矩阵,而矩阵对于我们而言,什么是矩阵?矩阵(Matrix)是一个按照长方阵列排列的复数或实数集合,最早来自于方程组的系数及常数所构成的方阵。这... -
cuda矩阵之心得
2017-07-18 09:35:45一、先了解什么是矩阵 在数学中,矩阵(Matrix)是一个按照长方阵列排列的复数或实数集合[1],最早来自于方程组的系数及常数所构成的方阵。 二、矩阵的基本参数 typedef struct { int width; int height; ... -
矩阵相乘优化算法实现讲解
2016-04-12 10:41:47矩阵是高等代数学中的常见工具,也常见于统计分析等应用数学学科中。并且在ACM竞赛,有很多涉及到矩阵知识的题。许多算法都会结合矩阵来处理,而比较具有代表性的矩阵算法有:矩阵快速幂、高斯消元等等。 -
可逆矩阵性质总结_线性代数下的行列式和矩阵
2021-01-01 23:36:10于是我们考虑的问题是:齐次方程组:是否存在非零解,以及存在的条件通解的结构与性质解法非齐次方程组:是否有解,以及有解的条件是什么有多少解以及对应解数量的条件是什么多解的结构与性质解法行列式二,三阶行列... -
机器学习笔记05-矩阵基础知识
2020-12-16 19:03:38在数学中,矩阵(Matrix)是一个按照长方阵列排列的复数或实数集合 ,最早来自于方程组的系数及常数所构成的方阵。这一概念由19世纪英国数学家凯利首先提出。直观解释就是由中括号包起来的一组二维数据,如下所示即... -
矩阵乘法实现(c语言版)
2016-05-09 12:41:58矩阵是高等代数学中的常见工具,也常见于统计分析等应用数学学科中。在物理学中,矩阵于电路学、力学、光学和量子物理中都有应用;计算机科学中,三维动画制作也需要用到矩阵。 矩阵的运算是数值分析领域的重要问题... -
洛谷 P1939 【模板】矩阵加速(数列):优化递推式的方法——矩阵快速幂
2017-12-15 19:49:08在大多数情况下,O(n)的效率都是值得骄傲的,然而,有时候并不是,比如如何...如果你还不知道矩阵快速幂是什么,请走这边:传送门 对于这道题,嗯,模板嘛,已经告诉你了式子,就只需要考虑矩阵了,对于整个过程,我们 -
数据结构和算法--Java实现矩阵
2018-12-09 15:10:43矩阵是高等代数学中的常见工具,也常见于统计分析等应用数学学科中。下面我们简单复习下。 什么是矩阵 1.矩阵定义 在数学中,矩阵(Matrix)是一个按照长方阵列排列的实数或复数的集合,最早来自于方程组的系数... -
矩阵初步(线性代数在信息学中的体现)
2019-08-06 17:20:35在数学上,矩阵是指纵横排列的二维数据表格,最早来自于方程组的系数及常数所构成的方阵。这一概念由19世纪英国数学家凯利首先提出。矩阵是高等代数学中的常见工具,也常见于统计分析等应用数学学科中。在物理学中,... -
睡前1小时数学之矩阵及其应用
2016-10-01 13:58:00线性方程组中的未知数的量排成一个矩阵,加上常数项,就是增广矩阵。 还有什么表示线性转换,就不讲了。反正。也不会。 2,标记。 一般将一个矩阵中的m*n个元素,简称元,数a(ij)的位于矩阵A的第i行,第j列。这... -
C++矩阵及其加速—————求斐波拉契数列第n项讲解
2019-04-02 13:33:38也许你只是不小心点入了此博客,为了不眛自己的良心,首先我们会介绍什么是矩阵。 概念: 在数学中,矩阵(Matrix)是一个按照长方阵列排列的实数或复数集合,最早来自于方程组的系数及常数所构成的方阵。 由m×n个... -
bzoj2738 矩阵乘法
2019-10-02 17:09:47【题解】 整体二分,然后用二维树状数组(单点修改区间查询)统计即可。...不知道为什么跑的特别慢(可能是我的整体二分常数太大?) # include <stdio.h> # include <string.h> # include &l... -
POJ -3233,HDU-5015(矩阵快速幂)
2019-03-29 21:50:23题目链接: POJ-3233 ...这次的转移矩阵是由一些子矩阵构成,不再是一些普通的常数了。 #include<cstdio> #include <iostream> #include <algorithm> #include <cstring> ... -
单应性矩阵和仿射变换_图像单应性变换理解
2021-01-13 12:46:15什么是单应性?图像中的2D点(x,y)可以被表示成3D向量的形式(x1,x2,x3),...单应变换矩阵是一个3*3的矩阵H。这个变换可以被任意乘上一个非零常数,而不改变变换本身。所以它虽然具有9个元素,但是具有8个自由度。这意... -
从特征值,特征向量,对角矩阵,特征值分解到奇异值分解
2020-12-07 10:37:59特征值,特征向量:设A是n阶矩阵,若存在常数λ及非零常数α,使得Aα=λα,称λ为矩阵A的一个特征值,α是属于特征值λ的矩阵A的一个特征向量。 (这样读起来好像没啥用,下面会给一个例子来说明怎么算 -
特征值_抽象矩阵的特征值与特征向量
2021-01-15 17:22:12//特征值与特征向量(4)//01前言(1)今天我们来讨论抽象矩阵的特征值与特征向量。这类抽象问题, 最重要的方法就是定义法!当然这里还需要一定的“观察力”, 这也是...注意f(x)若含有常数项a0, 则f(A)对应的项是a0E。... -
矩阵论(补充知识):特征多项式的展开式
2019-11-16 22:56:53国内线性代数教材上关于n阶矩阵AAA的特征多项式的系数只讲了常数项、n-1次项和n次项的,分别为det(A),tr(A),1det(A),tr(A),1det(A),tr(A),1。一直很好奇其他项的系数是什么样的。查资料知有如下定理: 定理:设A∈... -
特征值与特征向量_抽象矩阵的特征值与特征向量
2020-12-05 11:22:49//特征值与特征向量(4)//01前言(1)今天我们来讨论抽象矩阵的特征值与特征向量。这类抽象问题, 最重要的方法就是定义法!当然这里还需要一定的“观察力”, 这也是...注意f(x)若含有常数项a0, 则f(A)对应的项是a0E。... -
letcode-《二维区域和检索 - 矩阵不可变》
2021-03-03 00:27:40问题 问题:给定一个二维数组,给定一个矩形框的...数组已经是常数访问时间了,还有什么花哨的玩法么 代码 class NumMatrix { private int[][] matrix; public NumMatrix(int[][] matrix) { this.matrix = matrix; -
增广矩阵的秩判断解的个数_「线性代数」非齐次线性方程组,无解、唯一解、无穷解详细做法...
2020-12-31 08:24:14在线性代数中,我们往往会做到一类题目,那就是给定两个矩阵A、B,其中设有未知数,问我们什么时候AX=B无解、有唯一解、有无穷多解。我们认知中的AX=B便是非齐次线性方程组的表达式(常数项不全为零的线性方程组称为...