前言

回顾一下BLAS Level 1 2 3中的运算都有什么类型

• BLAS Level 1

  在BLAS Level 1中，进行的是向量-向量的操作。其中相关概念有

  • 向量类型: 实数域，复数域，共轭
  • 运算操作: 元素求和，向量相加，向量拷贝，点乘，二范数(欧几里得范数)，向量绕点旋转，Givens旋转，改进Givens旋转，计算Givens旋转参数，标量与向量的乘积，向量交换，查找最大最小值，复数域向量的绝对值

• BLAS Level 2

  在BLAS Level 2中，进行的是矩阵-向量的操作。其中相关概念有

  • 矩阵类型: 一般带状矩阵，一般矩阵，Hermitian带状矩阵，Hermitian矩阵，Hermitian压缩矩阵，对称带状矩阵，对称压缩矩阵，对称矩阵，三角带状矩阵，三角压缩矩阵，三角矩阵
  • 运算操作: 矩阵-向量乘积，一阶更新，二阶更新，解线性方程组

• BLAS Level 3

  在BLAS Level 3中，进行的是矩阵-矩阵的操作。其中相关概念有

  • 矩阵类型:一般矩阵，Hermitian矩阵，对称矩阵，三角矩阵
  • 运算操作: 矩阵乘法，k阶更新，2k阶更新， 解方程

向量矩阵类型学习

向量

• 实数域
  

  x=[1 , 2 , 3.1 , 4]
  

  x=[1\\ ,\\ 2\\ ,\\ 3.1\\ ,\\ 4]
  
  元素为实数的就叫实数域上的向量

• 复数域
  

  x=[−i,−2i,−3i,−4i]x=[4−i , 3−2i , 2−3i , 1−4i]
  

  x=[-i , -2i,-3i,-4i]\\\\
  x=[4-i \\ ,\\ 3-2i\\ ,\\ 2-3i\\ ,\\ 1-4i]
  
  元素为复数的就叫复数域上的向量，第一个向量的实部为0，第二个向量的实部非0

• 共轭向量

  当两个复数的实部相等，虚部互为相反数时，这两个复数叫做互为共轭复数，实数的共轭复数就是本身

矩阵

• 一般矩阵
  

  A=⎡⎣⎢a11a21a31a12a22a32a13a23a33⎤⎦⎥
  

  A=\\begin{bmatrix}
  a_{11}&a_{12}&a_{13}\\\\
  a_{21}&a_{22}&a_{23}\\\\
  a_{31}&a_{32}&a_{33}
  \\end{bmatrix}

• 带状矩阵

  如果矩阵

  A=aij
  A=a_{ij}是一个

  n∗n
  n*n的矩阵，假设在对角边界带状区域之外的所有元素都是0，其中边界范围由常数

  k1
  k_1和

  k2
  k_2决定，那么就有
  

  当j<(i−k1)或者j>(i+k2)的时候有aij=0
  

  当j(i+k_2)的时候有a_{ij}=0
  
  这里面的

  k1
  k_1称为下带宽(lower bandwidth)，

  k2
  k_2称为上带宽(upper bandwidth)
  

  A=⎡⎣⎢a1100a12a2200a23a33⎤⎦⎥
  

  A=\\begin{bmatrix}
  a_{11}&a_{12}&0\\\\
  0&a_{22}&a_{23}\\\\
  0&0&a_{33}
  \\end{bmatrix}

• 三角带状矩阵

  当带状矩阵的

  k1=k2=1
  k_1=k_2=1的时候，就是三角带状矩阵了
  

  A=⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢B11B210⋮⋮0B12B22B32⋱⋱⋯0B23B33B43⋱⋯⋯⋱B34B44B540⋯⋱⋱B45B55B650⋮⋮0B56B66⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥
  

  A=\\begin{bmatrix}
  B_{11}&B_{12}&0&\\cdots&\\cdots&0\\\\
  B_{21}&B_{22}&B_{23}&\\ddots&\\ddots&\\vdots\\\\
  0&B_{32}&B_{33}&B_{34}&\\ddots&\\vdots\\\\
  \\vdots&\\ddots&B_{43}&B_{44}&B_{45}&0\\\\
  \\vdots&\\ddots&\\ddots&B_{54}&B_{55}&B_{56}\\\\
  0&\\cdots&\\cdots&0&B_{65}&B_{66}
  \\end{bmatrix}
  
  ​

• Hermitian矩阵

  Hermitian矩阵(也称为self-adjoint矩阵)，是一个复数方阵，与它自己的共轭转置是相等的；也就是说第i行第j列的元素等于它的共轭转置矩阵的第j行第i列，对于所有的索引i,j都有
  

  aij=aji¯¯¯¯A=AH¯¯¯¯¯
  

  a_{ij}=\\overline{a_{ji}}\\\\
  A=\\overline{A^H}
  
  一个实例就是
  

  A=⎡⎣⎢22−i42+i3−i4i1⎤⎦⎥
  

  A=\\begin{bmatrix}
  2&2+i&4\\\\
  2-i&3&i\\\\
  4&-i&1
  \\end{bmatrix}

• 三角矩阵
  

  A=⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢a11a21a31⋮an1a22a32⋮an2⋱⋱⋯⋱an,n−10an,n⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥
  

  A=\\begin{bmatrix}
  a_{11}& & & & 0\\\\
  a_{21}&a_{22}& & &\\\\
  a_{31}&a_{32}&\\ddots& &\\\\
  \\vdots&\\vdots&\\ddots&\\ddots& \\\\
  a_{n1}&a_{n2}&\\cdots&a_{n,n-1}&a_{n,n}
  \\end{bmatrix}
  
  称为下三角矩阵(lower triangular matrix)或者左三角矩阵(left triangular matrix)。同理还有上三角矩阵就不说了。

• 对称矩阵

  对称矩阵就是一个转置等于它自身的一个方阵
  

  A=AT
  

  A=A^T
  
  例如
  

  A=⎡⎣⎢17374−53−56⎤⎦⎥
  

  A=\\begin{bmatrix}
  1&7&3\\\\
  7&4&-5\\\\
  3&-5&6&
  \\end{bmatrix}

• 压缩存储

  在Cblas中提供了很多压缩存储，但是目前还没分析透彻，有兴趣可以看看英文文档，其中提供了各种转换压缩存储的主程序。主要包含的压缩存储有:

  • 带状压缩存储
  • 按列存储时候Hermitian矩阵的上三角带状存储；按行存储时候Hermitian矩阵的下三角带状存储
  • 对称矩阵的上三角压缩存储；对称矩阵的下三角压缩存储
  • 按列存储时候上三角带状矩阵或者下三角带状矩阵的压缩存储；按行存储的时候上三角带状矩阵或者下三角带状矩阵的压缩存储

向量矩阵操作学习

向量-向量

• Givens旋转

  直接从维基百科上搬过来：

  == 矩阵表示 ==
  吉文斯旋转表示为如下形式的[[矩阵]]
  

  G(i,j,θ)=⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢1⋮0⋮0⋮0⋯⋱⋯⋯⋯0⋮c⋮s⋮0⋯⋯⋱⋯⋯0⋮−s⋮c⋮0⋯⋯⋯⋱⋯0⋮0⋮0⋮1⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥
  

  G(i, j, \\theta) = \\begin{bmatrix} 1 & \\cdots & 0 & \\cdots & 0 & \\cdots & 0 \\\\
  \\vdots & \\ddots & \\vdots & & \\vdots & & \\vdots \\\\
  0 & \\cdots & c & \\cdots & -s & \\cdots & 0 \\\\
  \\vdots & & \\vdots & \\ddots & \\vdots & & \\vdots \\\\
  0 & \\cdots & s & \\cdots & c & \\cdots & 0 \\\\
  \\vdots & & \\vdots & & \\vdots & \\ddots & \\vdots \\\\
  0 & \\cdots & 0 & \\cdots & 0 & \\cdots & 1
  \\end{bmatrix}
  
  这里的 c= cos(θ) 和 s = sin(θ) 出现在第 ”i” 行和第 ”j” 行与第 ”i” 列和第 ”j” 列的交叉点上。就是说，吉文斯旋转矩阵的所有非零元定义如下：
  

  gkkgiigjjgijgji=1fork≠i,j=c=c=−s=s
  

  \\begin{align}
  g_{k\\, k} &{}= 1 \\qquad \\text{for} k \\neq i,j \\\\
  g_{i\\, i} &{}= c \\\\
  g_{j\\, j} &{}= c \\\\
  g_{i\\, j} &{}= -s \\\\
  g_{j\\, i} &{}= s
  \\end{align}

  乘积

  G(i, j, θ)∗x
  G(i, j, θ)*x 表示向量

  x
  x在 (i,j)(i,j)平面中的逆时针旋转 θ 弧度。

  Givens 旋转在数值线性代数中主要的用途是在向量或矩阵中介入零。例如，这种效果可用于计算矩阵的 QR分解。超过Householder变换的一个好处是它们可以轻易的并行化，另一个好处是对于非常稀疏的矩阵计算量更小。

  == 稳定计算 ==
  当一个吉文斯旋转矩阵

  G(i,j,θ)
  G(i,j,θ)从左侧乘另一个矩阵

  A
  A的时候，GAGA 只作用于

  A
  A的第 ii和

  j
  j 行。所以我们集中于下列问题。给出 aa和

  b
  b，找到 c=cos(θ)c = cos(θ) 和

  s=sin(θ)
  s= sin (θ) 使得
  

  [cs−sc][ab]=[r0]
  

  \\begin{bmatrix}
  c & -s \\\\
  s & c \\end{bmatrix}
  \\begin{bmatrix}
  a \\\\ b
  \\end{bmatrix}=
  \\begin{bmatrix}
  r \\\\ 0
  \\end{bmatrix}
  
  明确计算出 θ 是没有必要的。我们转而直接获取 ”c”, ”s” 和 ”r”。一个明显的解是
  

  rcs=a2+b2−−−−−−√=a/r=−b/r
  

  \\begin{align}
  r &= \\sqrt{a^2 + b^2} \\\\
  c &= a / r \\\\
  s &=-b / r
  \\end{align}

矩阵-向量

• Hermitian矩阵或者其压缩存储形式

  • 一阶更新

  A:=α∗x∗conjg(x′)+A
  

  A:=\\alpha*x*conjg(x\’)+A

  其中

  x
  x是一个nn维向量

  • 二阶更新

  A:=α∗x∗conjg(y′)+conjg(α)∗y∗conjg(x′)+A
  

  A := \\alpha *x*conjg(y\’) + conjg(\\alpha)*y *conjg(x\’) + A

  其中

  x,y
  x,y都是

  n
  n向量

• 对称矩阵或者其压缩存储形式

  • 一阶更新
  a:=α∗x∗x′+A

  a:= \\alpha*x*x\’+ A

  ​ 其中

  x
  x为nn维矩阵

  • 二阶更新
    

    A:=α∗x∗y′+α∗y∗x′+A
    

    A:= \\alpha*x*y\’+ \\alpha*y*x\’ + A
    
    其中

    x,y
    x,y都是

    n
    n维矩阵

矩阵-矩阵

• Hermitian矩阵

  • k阶更新
    

    C:=α∗A∗AH+β∗C

    C := \\alpha*A*A^H + \\beta*C
    
    其中

    A
    A是n∗kn*k维矩阵，

    C
    C是n∗nn*n的Hermitian矩阵

  • 2k阶更新
    

    C:=α∗A∗BH+conjg(α)∗B∗AH+β∗C,
    

    C := \\alpha*A*B^H + conjg(\\alpha)*B*A^H + \\beta*C,
    
    其中

    A,B
    A,B都是

    n∗k
    n*k的矩阵，

    C
    C是n∗nn*n的Hermitian矩阵

• 对称矩阵

  • k阶更新

  C:=α∗A∗A′+β∗C
  

  C := \\alpha*A*A\’ + \\beta*C

  ​ 其中

  A
  A是n∗kn*k的矩阵，

  C
  C是n∗nn*n的矩阵

  • 2k阶更新
    

    C:=α∗A∗B′+α∗B∗A′+β∗C
    

    C := \\alpha*A*B\’ + \\alpha*B*A\’ + \\beta*C
    
    其中

    A,B
    A,B都是

    n∗k
    n*k的矩阵，

    C
    C是n∗nn*n的矩阵

后续

目前感觉也就这么多内容了。关于后面的稀疏矩阵计算相关内容，等用到了再学习。后面着重实现以下最基本的全矩阵(full matrix)运算，也就是说不考虑矩阵是对称的、Hermitian的或者带状什么的。主要包含向量-向量、矩阵-向量、矩阵-矩阵的乘法关系。