2.2 矩阵
矩阵可以用来表示数据集，描述数据集上的变换，是MADlib中数据的基本格式，通常使用二维数组数据类型。MADlib中的向量是一维数组，可看作是矩阵的一种特殊形式。MADlib的矩阵运算模块（matrix_ops）实现SQL中的矩阵操作。本节将介绍矩阵的概念，说明MADlib矩阵运算相关函数，并举出一些简单的函数调用示例。
2.2.1 矩阵定义
矩阵（matrix）是把数集合汇聚成行和列的一种表表示。术语“m×n矩阵”通常用来说明矩阵具有m行和n列。下面所示的矩阵A是2×3矩阵。如果m=n，则我们称该矩阵为方阵（square matrix）。矩阵A的转置记作AT（通过交换A的行和列得到）。
矩阵的元素用带下标的小写字母表示。对于矩阵A，aij是其第i行第j列的元素。行自上而下编号，列自左向右编号，编号从1开始。例如，a21是矩阵A的第2行第1列的元素，该元素的值是7。
矩阵的每一行或列定义一个向量。对于矩阵A，其第i个行向量（row vector）可以用ai∗表示，第j个列向量（column vector）用a∗j表示。使用前面的例子，a2∗=[7 5 2]，而a∗3=[1 2]T。注意，行向量和列向量都是矩阵，必须加以区分，即元素个数相同并且值相同的行向量和列向量代表不同的矩阵。
2.2.2 MADlib中的矩阵表示
MADlib支持稠密和稀疏两种矩阵表示形式，所有矩阵运算都以任一种表示形式工作。
1. 稠密
矩阵被表示为一维数组的行集合，例如3×10的矩阵如下：
row_id列表示每一行的行号，是从1到N没有重复值的连续整型序列，N为矩阵的行数。row_vec列对应构成矩阵每行的一维数组，即行向量。
2. 稀疏
使用行列下标指示矩阵中每一个非零项，例如：
常用这种方式表示包含多个零元素的稀疏矩阵。上面的例子只用6行表示一个4×7的矩阵中的非零元素。矩阵的行列元素个数分别由row_id和col_id的最大值指定。注意最后一行，即使value为0也要包含此行，它指出了矩阵的维度，而且指示矩阵的第4行与第7列的元素值都是0。
对于稀疏矩阵表，row_id和col_id列逻辑类似于关系数据库的联合主键，要求非空且唯一。value列应该是标量简单数据类型。上面矩阵对应的稠密表示如下：
2.2.3 MADlib中的矩阵运算函数
与向量操作相同，矩阵运算函数支持的元素数据类型也包括SMALLINT、INTEGER、BIGINT、FLOAT8和NUMERIC（内部被转化为FLOAT8，可能丢失精度）。
1. 矩阵操作函数分类
MADlib的矩阵操作函数可分为表示、计算、提取、归约、创建、转换六类。下面列出每一类中所包含的函数名称及其参数。
（1）表示函数
-- 转为稀疏矩阵 matrix_sparsify( matrix_in, in_args, matrix_out, out_args) -- 转为稠密矩阵 matrix_densify( matrix_in, in_args, matrix_out, out_args) -- 获取矩阵的维度 matrix_ndims( matrix_in, in_args )
（2）计算函数
-- 矩阵转置 matrix_trans( matrix_in, in_args, matrix_out, out_args) -- 矩阵相加 matrix_add( matrix_a, a_args, matrix_b, b_args, matrix_out, out_args) -- 矩阵相减 matrix_sub( matrix_a, a_args, matrix_b, b_args, matrix_out, out_args) -- 矩阵相乘 matrix_mult( matrix_a, a_args, matrix_b, b_args, matrix_out, out_args) -- 数组元素相乘 matrix_elem_mult( matrix_a, a_args, matrix_b, b_args, matrix_out, out_args) -- 标量乘矩阵 matrix_scalar_mult( matrix_in, in_args, scalar, matrix_out, out_args) -- 向量乘矩阵 matrix_vec_mult( matrix_in, in_args, vector)
（3）提取函数
-- 从行下标提取行 matrix_extract_row( matrix_in, in_args, index) -- 从列下标提取列 matrix_extract_col( matrix_in, in_args, index) -- 提取主对角线元素 matrix_extract_diag( matrix_in, in_args)
（4）归约函数（指定维度的聚合）
-- 获取指定维度的最大值。如果fetch_index=True，返回对应的下标 matrix_max( matrix_in, in_args, dim, matrix_out, fetch_index) -- 获取指定维度的最小值。如果fetch_index=True，返回对应的下标 matrix_min( matrix_in, in_args, dim, matrix_out, fetch_index) -- 获取指定维度的和 matrix_sum( matrix_in, in_args, dim) -- 获取指定维度的均值 matrix_mean( matrix_in, in_args, dim) -- 获取矩阵范数 matrix_norm( matrix_in, in_args, norm_type)
（5）创建函数
-- 创建一个指定行列维度的矩阵，用1初始化元素值 matrix_ones( row_dim, col_dim, matrix_out, out_args) -- 创建一个指定行列维度的矩阵，用0初始化元素值 matrix_zeros( row_dim, col_dim, matrix_out, out_args) -- 创建单位矩阵 matrix_identity( dim, matrix_out, out_args) -- 用给定对角元素初始化矩阵 matrix_diag( diag_elements, matrix_out, out_args)
（6）转换函数
-- 矩阵求逆 matrix_inverse( matrix_in, in_args, matrix_out, out_args) -- 广义逆矩阵 matrix_pinv( matrix_in, in_args, matrix_out, out_args) -- 矩阵特征提取 matrix_eigen( matrix_in, in_args, matrix_out, out_args) -- Cholesky分解 matrix_cholesky( matrix_in, in_args, matrix_out_prefix, out_args) -- QR分解 matrix_qr( matrix_in, in_args, matrix_out_prefix, out_args) -- LU分解 matrix_lu( matrix_in, in_args, matrix_out_prefix, out_args) -- 求矩阵的核范数 matrix_nuclear_norm( matrix_in, in_args) -- 求矩阵的秩 matrix_rank( matrix_in, in_args)
矩阵转换函数仅基于内存操作实现。单一节点的矩阵数据被用于分解计算。这种操作只适合小型矩阵，因为计算不是分布到多个节点执行的。
2. 矩阵操作函数示例
执行下面的脚本创建两个稠密表示的矩阵测试表并添加数据。mat_a矩阵4行4列，mat_b矩阵5行4列。
drop table if exists mat_a;
create table mat_a (row_id integer, row_vec integer[]);
insert into mat_a (row_id, row_vec) values
(1, '{9,6,5,8}'), (2, '{8,2,2,6}'), (3, '{3,9,9,9}'), (4, '{6,4,2,2}');
drop table if exists mat_b;
create table mat_b (row_id integer, vector integer[]);
insert into mat_b (row_id, vector) values
(1, '{9,10,2,4}'), (2, '{5,3,5,2}'), (3, '{0,1,2,3}'), (4, '{2,9,0,4}'), (5,
'{3,8,7,7}');
（1）由稠密矩阵表生成稀疏表示的表
drop table if exists mat_a_sparse;
select madlib.matrix_sparsify('mat_a', 'row=row_id, val=row_vec',
'mat_a_sparse', 'col=col_id, val=val');
drop table if exists mat_b_sparse;
select madlib.matrix_sparsify('mat_b', 'row=row_id, val=vector',
'mat_b_sparse', 'col=col_id, val=val');
madlib.matrix_sparsify函数将稠密表示矩阵表转为稀疏表示的矩阵表，四个参数分别指定输入表名、输入表参数（代表行ID的列名、存储矩阵元素值的列名等）、输出表名、输出表参数（代表列ID的列名、存储矩阵元素值的列名等）。
上面的例子将稠密矩阵转为稀疏表示，并新建表存储转换结果。源表的两列类型分别是整型和整型数组，输出表包含三列，行ID列名与源表相同，列ID列和值列由参数指定。由于mat_a表的矩阵中不存在0值元素，生成的稀疏矩阵表共有16条记录，而mat_b中有两个0值，因此稀疏表中只有18条记录。
（2）矩阵转置
matrix_trans函数的第一个参数是源表名，第二个参数指定行、列或值的字段名，第三个参数为输出表名。
-- 稠密格式
drop table if exists mat_a_r;
select madlib.matrix_trans('mat_a', 'row=row_id, val=row_vec','mat_a_r');
select * from mat_a_r order by row_id;
结果：
-- 稀疏格式
drop table if exists mat_b_sparse_r;
select madlib.matrix_trans('mat_b_sparse', 'row=row_id, col=col_id,
val=val','mat_b_sparse_r');
select * from mat_b_sparse_r order by row_id, col_id;
结果：
源矩阵5行4列，转置后的矩阵为4行5列。
（3）提取矩阵的主对角线
select madlib.matrix_extract_diag('mat_b', 'row=row_id, val=vector'),
madlib.matrix_extract_diag
('mat_b_sparse_r', 'row=row_id, col=col_id, val=val');
结果：
matrix_extract_diag函数的返回值是由对角线元素组成的数组。可以看到，矩阵和其对应的转置矩阵具有相同的主对角线。也就是说，矩阵转置实际上是沿着主对角线的元素对折操作。
（4）创建对角矩阵
结果：
madlib.matrix_diag函数输出的是一个稀疏表示的对角矩阵表，如果不指定“col=col_id”，输出表中代表列的列名为col。
（5）创建单位矩阵
drop table if exists mat_r; select madlib.matrix_identity(4, 'mat_r'); select * from mat_r;
结果：
matrix_identity函数创建一个稀疏表示的单位矩阵表。主对角线上的元素都为1、其余元素全为0的方阵称为单位矩阵。
（6）提取指定下标的行或列
结果返回两个向量，即mat_a的第2行、mat_b_sparse的第3列：
（7）获取指定维度的最大最小值及其对应的下标
drop table if exists mat_max_r, mat_min_r;
select madlib.matrix_max
('mat_a', 'row=row_id, val=row_vec', 2, 'mat_max_r', true),
madlib.matrix_min
('mat_b_sparse', 'row=row_id, col=col_id', 1, 'mat_min_r', true);
select * from mat_max_r, mat_min_r;
结果：
matrix_max和matrix_min函数分别返回指定维度的最大值和最小值，其中维度参数的取值只能是1或2，分别代表行和列。返回值为数组类型，如果最后一个参数为‘true’，表示结果中包含最大最小值对应的下标数组列。
（8）创建元素为全0的矩阵
drop table if exists mat_r01, mat_r02;
select madlib.matrix_zeros(3, 2, 'mat_r01', 'row=row_id, col=col_id,
val=entry'),
madlib.matrix_zeros(3, 2, 'mat_r02', 'fmt=dense');
select * from mat_r01;
select * from mat_r02;
结果分别为：
注意，元素值全为0，所以稀疏表示的矩阵表只有1行。
（9）创建元素为全1的矩阵
drop table if exists mat_r11, mat_r12;
select madlib.matrix_ones(3, 2, 'mat_r11', 'row=row_id, col=col_id,
val=entry'),
madlib.matrix_ones(3, 2, 'mat_r12', 'fmt=dense');
select * from mat_r11 order by row_id;
select * from mat_r12 order by row;
结果分别为：
因为元素值全为1，所以稀疏表示的矩阵表有6行。
（10）获取行列维度数
select madlib.matrix_ndims('mat_a', 'row=row_id, val=row_vec'),
madlib.matrix_ndims('mat_a_sparse', 'row=row_id, col=col_id');
结果：
（11）矩阵相加
与向量一样，矩阵也可以通过将对应元素（分量）相加来求和。MADlib的矩阵相加函数要求两个矩阵具有相同的行数和列数。更明确地说，假定A和B都是m×n的矩阵，A和B的和是m×n矩阵C，其元素由下式计算：
cij=aij+bij
结果：
madlib.matrix_add函数有三组参数，分别是两个相加的矩阵表和结果矩阵表。相加的两个矩阵表不必有相同的表示形式，如上面的函数调用中，一个矩阵为稠密形式，一个矩阵为稀疏形式，但两个矩阵必须具有相同的行列数，否则会报如下错误：
Matrix error: The dimensions of the two matrices don't match
矩阵加法具有如下性质：
•矩阵加法的交换律。加的次序不影响结果：A+B=B+A。
•矩阵加法的结合律。相加时矩阵分组不影响结果：(A+B)+C=A+(B+C)。
•矩阵加法单位元的存在性。存在一个零矩阵（zero matrix），其元素均为0并简记为0，是单位元。对于任意矩阵A，有A+0=A。
•矩阵加法逆元的存在性。对于每个矩阵A，都存在一个矩阵-A，使得A+(-A)=0。-A的元素为-aij。
（12）标量与矩阵相乘
与向量一样，也可以用标量乘以矩阵。标量α和矩阵A的乘积是矩阵B=αA，其元素由下式给出：
bij=αaij
例如，下面matrix_scalar_mult函数执行结果是由原矩阵的每个元素乘以3构成的矩阵表。
drop table if exists mat_r;
select madlib.matrix_scalar_mult('mat_a', 'row=row_id, val=row_vec', 3,
'mat_r');
select * from mat_r order by row_id;
结果：
矩阵的标量乘法具有与向量的标量乘法非常相似的性质。
•标量乘法的结合律，被两个标量乘的次序不影响结果：α(βA)=(αβ)A。
•标量加法对标量与矩阵乘法的分配率，两个标量相加后乘以一个矩阵等于每个标量乘以该矩阵之后的结果矩阵相加：(α+β)A=αA+βA。
•标量乘法对矩阵加法的分配率，两个矩阵相加之后的和与一个标量相乘等于每个矩阵与该标量相乘然后相加：α(A+B)=αA+αB。
•标量单位元的存在性，如果α=1，则对于任意矩阵A，有αA=A。
我们可以认为矩阵由行向量或列向量组成，因此矩阵相加或用标量乘以矩阵等于对应行向量或列向量相加或用标量乘它们。
（13）矩阵乘法
我们可以定义矩阵的乘法运算。先定义矩阵与向量的乘法。矩阵与列向量的乘法：m×n矩阵A乘以n×1的列矩阵u的积是m×1的列矩阵v=Au，其元素由下式给出：
vi=ai∗·uT
换言之，我们取A的每个行向量与u的转置的点积。注意，在下面的例子中，u的行数必然与A的列数相等。
类似地，我们可以定义矩阵被行向量左乘。矩阵与行向量的乘法：1×m的行矩阵u乘以m×n矩阵A的积是1×n的行矩阵v=uA，其元素由下式给出：
vi=u·(a∗j)T
换言之，我们取该行向量与矩阵A的每个列向量的转置的点积。下面给出一个例子：
MADlib的matrix_vec_mult函数用于计算一个m×n矩阵乘以一个1×n的矩阵（向量），结果是一个1×m的矩阵。如下面的5×4的矩阵mat_b乘以一个1×4的矩阵，结果是一个1×5的矩阵。
可以用下面的查询验证矩阵乘以向量的结果。
dm=# select array_agg(madlib.array_dot(vector,array[1,2,3,4])) from mat_b;
array_agg
------------------
{51,34,20,36,68}
(1 row)
我们定义两个矩阵的乘积，作为上述概念的推广。m×n矩阵A与n×p矩阵B的积是m×p矩阵C（C=AB），其元素由下式给出：
cij=ai∗·(b∗j)T
换言之，C的第ij个元素是A的第i个行向量与B的第j个列向量转置的点积。
matrix_mult函数用于矩阵相乘。如前所述，第一组参数中的矩阵列数应该与第二组参数中的矩阵行数相同，否则会报错：
可以对mat_b先进行转置，再与mat_a相乘。matrix_mult函数调用时的trans=true参数表示先对mat_b表行列转置再进行矩阵乘法。这次的矩阵乘法计算将正常执行。
结果是一个4×5矩阵：
执行结果与下面的查询相同。
矩阵乘法具有如下性质。
•矩阵乘法的结合律，矩阵乘的次序不影响计算结果：(AB)C=A(BC)。
•矩阵乘法的分配率，矩阵乘法对矩阵加法是可分配的：A(B+C)=AB+AC并且(B+C)A=BA+CA。
•矩阵乘法单位元的存在性，若Ip是p×p矩阵的单位矩阵，则对于任意m×n矩阵A，AIn=A并且ImA=A。
一般地，矩阵乘法是不可交换的，即AB≠BA。
如果我们有一个n×1列向量u，我们就可以把m×n矩阵A被该向量右乘看作u到m维列向量v=Au的变换。类似地，如果我们用一个（行）向量u=[u1,…,um]左乘A，我们可以将它看作u到n维行向量v=uA的变换。这样，我们可以把一个任意m×n矩阵A看作一个把一个向量映射到另一个向量空间的函数。
（14）两矩阵元素相乘
与矩阵乘法定义不同，MADlib的两矩阵元素相乘定义为C=AB，A、B、C均为m×n矩阵，C的元素由下式给出：
cij=aij×bij
MADlib的matrix_elem_mult函数执行两矩阵元素相乘，并输出结果矩阵。
结果：
（15）求矩阵的秩
select madlib.matrix_rank('mat_a', 'row=row_id, val=row_vec');
结果：
注意，当矩阵以稀疏形式表示，并且列数大于行数时，matrix_rank函数会报错。
dm=# select madlib.matrix_rank('mat_b_sparse_r', 'row=row_id, col=col_id,
val=val');
ERROR: plpy.SPIError: Function
"madlib.__matrix_compose_sparse_transition(double
precision[],integer,integer,integer,integer,double precision)": Invalid col
id. (UDF_impl.hpp:210) (seg20 hdp4:40000 pid=123035) (plpython.c:4663)
CONTEXT: Traceback (most recent call last):
PL/Python function "matrix_rank", line 23, in <module>
return matrix_ops.matrix_rank(schema_madlib, matrix_in, in_args)
PL/Python function "matrix_rank", line 2702, in matrix_rank
PL/Python function "matrix_rank", line 2672, in matrix_eval_helper
PL/Python function "matrix_rank"
dm=#
矩阵的秩（rank of a matrix）常常用来刻画矩阵。设矩阵A=(aij)m×n，在A中任取k行k列交叉处元素按原相对位置组成的k阶行列式，称为A的一个k阶子式。m×n矩阵A共有
个k阶子式。若A有r阶子式不为0，任何r+1阶子式（如果存在的话）全为0，则称r为矩阵A的秩，记作R(A)。
矩阵的秩具有以下基本性质：
•0矩阵的秩为0。
•若R(A)=r，则A中至少有一个r阶子式Dr≠0，所有r+1阶子式为0，且更高阶子式均为0，r是A中非零子式的最高阶数。
•矩阵转置，秩不变。
•0≤R(A)≤min(m,n)。
•若A是n×n方阵，并且|A|≠0，则R(A)=n；反之，若R(A)=n，则|A|≠0。
矩阵的秩（rank of a matrix）是行空间和列空间的最小维度，此维度中的向量组是线性无关的。例如，把一个1×n的行向量复制m次，产生一个m×n的矩阵，则我们只有一个秩为1的矩阵。
（16）求逆矩阵
drop table if exists mat_r;
select madlib.matrix_inverse('mat_a', 'row=row_id, val=row_vec', 'mat_r');
select row_vec from mat_r order by row_id;
结果：
设A、B是两个矩阵，若AB=BA=E，则称B是A的逆矩阵，而A则被称为可逆矩阵。其中E是单位矩阵。下面看一个不可逆矩阵的例子。
create table t1 (a int, b int[]);
insert into t1 values
(1,'{1,2,3}'),(2,'{2,4,6}'),(3,'{3,6,9}');
select madlib.matrix_rank('t1', 'row=a, val=b');
select madlib.matrix_inverse('t1', 'row=a, val=b', 't2');
select * from t2 order by a;
3阶矩阵t1的秩为1，用matrix_inverse求t1的逆矩阵，结果如下：
如果求逆的矩阵不是方阵，那matrix_inverse函数会报如下错误：
Matrix error: Inverse operation is only defined for square matrices
（17）求广义逆矩阵
把逆矩阵推广到不可逆方阵（奇异矩阵）或长方矩阵上，这就是所谓的广义逆矩阵。广义逆矩阵具有逆矩阵的部分性质，并且在方阵可逆时，它通常与逆矩阵一致。
drop table if exists mat_r;
select madlib.matrix_pinv('mat_a', 'row=row_id, val=row_vec', 'mat_r');
select row_vec from mat_r order by row_id;
结果：
matrix_pinv函数用于求矩阵的广义逆矩阵。还以上面的不可逆方阵为例，求它的广义逆矩阵。
drop table if exists t1,t2;
create table t1 (a int, b int[]);
insert into t1 values
(1,'{1,2,3}'),(2,'{2,4,6}'),(3,'{3,6,9}');
select madlib.matrix_pinv('t1', 'row=a, val=b', 't2');
select * from t2 order by a;
结果：
再看一个长方矩阵的例子。
drop table if exists mat_r;
select madlib.matrix_ndims('mat_b', 'row=row_id, val=vector'),
madlib.matrix_pinv('mat_b', 'row=row_id, val=vector', 'mat_r');
select * from mat_r order by row_id;
mat_b是一个5×4矩阵，它的广义逆矩阵如下：
（18）提取矩阵的特征值
drop table if exists mat_r;
select madlib.matrix_eigen('mat_a', 'row=row_id, val=row_vec', 'mat_r');
select * from mat_r order by row_id;
结果：
（19）求矩阵范数
matrix_norm函数用于求矩阵范数，支持的类型值有‘fro’‘one’‘inf’‘max’‘spec’，分别代表frobenius范数、1范数、infinity范数、max范数和spectral范数。默认为frobenius范数。
select madlib.matrix_norm('mat_b_sparse', 'row=row_id, col=col_id, val=val');
结果：
F-范数的公式为：
。依据公式可知下面查询的结果与matrix_norm函数的返回值相等。
select sqrt(sum(power(val,2))) from mat_b_sparse;
（20）求矩阵核范数
select madlib.matrix_nuclear_norm('mat_a', 'row=row_id, val=row_vec');
结果：
矩阵的核范数是指矩阵奇异值的和，关于矩阵奇异值，在讨论MADlib的矩阵分解函数时再进行详细说明。