矩阵运算
在处理面板数据时,矩阵运算相对于表运算而言,速度更快,效率更高。DolphinDB 针对矩阵运算场景为用户提供了丰富且便捷的内置函数,其中部分矩阵运算采用 OpenBlas 和 Lapack 进行优化,与 Matlab 的性能在一个数量级。 本教程将重点介绍 DolphinDB 矩阵的特点及矩阵运算的相关方法。
本教程将讲述以下有关矩阵内容:
1. 矩阵的存储
DolphinDB 中的矩阵采用列优先(column major)的形式。因此创建矩阵时,数据在内存中按列存储。系统会将输入的向量数据从左向右以列为单位填充为矩阵中的元素。注意:矩阵中行和列的下标都是从0开始的。
例:用[1, 2, 3], [4, 5, 6]两个向量创建矩阵得到的矩阵维度为 3*2,而不是 2*3。
>matrix(1 2 3, 4 5 6);
#0 #1
-- --
1 4
2 5
3 6
由于 DolphinDB 的矩阵是列优先的,读取或构建矩阵时,列操作比行操作要更为高效。
给定一个长度为1000的向量 v,构建一个900*100的矩阵。其第 i 行是 v[i:i+100]。比较下列两种做法:
def alongRows(v, mutable m){
for(i in 0:900){
m[i,]=v[i:(i+100)]
}
}
v=1..1000
m = matrix(int,900,100)
timer(10){alongRows(v, m)}
上述脚本逐行对矩阵进行赋值,耗时为236.369毫秒。
def alongColumns(v, mutable m){
for(i in 0:100){
m[,i]=v[i:(i+900)]
}
}
v=1..1000
m = matrix(int,900,100)
timer(10){alongColumns(v, m)}
上述脚本逐列对矩阵进行赋值,耗时为2.991毫秒。对比按行赋值,耗时缩短了两个数量级。
2. 矩阵的基础操作
2.1. 矩阵的创建
DolphinDB 的矩阵有三种类型,分别为普通矩阵(matrix)、索引序列(indexed series)和索引矩阵(indexed matrix)。
普通矩阵:调用函数
matrix
创建一个指定数据类型和维度的普通矩阵。函数 matrix 也可根据已有的数据包括向量、矩阵、表、元组及这些数据结构的组合创建一个矩阵。索引序列:索引序列带行索引标签的单列矩阵,可以通过函数
indexedSeries
进行创建,或者通过一个单列的矩阵调用setIndexedSeries
函数进行转换。索引矩阵:索引矩阵为带行/列索引标签的矩阵。只能通过普通矩阵调用函数
setIndexedMatrix
生成。
除直接生成矩阵外,通过运算符 cast($)
可以把一个向量重组为一个矩阵。
>m=1..10$5:2;
>m;
#0 #1
-- --
1 6
2 7
3 8
4 9
5 10
函数 reshape
可以实现矩阵与向量之间的互相转换。
>m=(1..6).reshape(3:2);
>m;
#0 #1
-- --
1 4
2 5
3 6
>y=m.reshape()
>y;
[1,2,3,4,5,6]
reshape 实现的功能相当于 cast($)
和 flatten
的组合。函数 flatten
可把一个矩阵转换为一个向量。
>x = flatten m;
>x
[1,2,3,4,5,6]
>x$2:5
#0 #1 #2 #3 #4
-- -- -- -- --
1 3 5 7 9
2 4 6 8 10
此外,可以调用函数生成一些特殊的矩阵类型:
单位矩阵
使用函数
eye(n)
创建维度为 n 的单位矩阵。>eye(3); #0 #1 #2 -- -- -- 1 0 0 0 1 0 0 0 1
对角矩阵
函数
diag(X)
:如果 X 是一个向量,返回一个对角矩阵,X 为主对角线上的元素;如果 X 是一个方阵,返回一个由主对角线元素组成的向量。>m=diag(1..5); >m; #0 #1 #2 #3 #4 -- -- -- -- -- 1 0 0 0 0 0 2 0 0 0 0 0 3 0 0 0 0 0 4 0 0 0 0 0 5 >diag(m); [1,2,3,4,5]
使用rename!
函数给矩阵添加行标签和列标签:
>m=1..9$3:3;
>m;
#0 #1 #2
-- -- --
1 4 7
2 5 8
3 6 9
>m.rename!(`col1`col2`col3);
col1 col2 col3
---- ---- ----
1 4 7
2 5 8
3 6 9
>m.rename!(1 2 3, `c1`c2`c3);
c1 c2 c3
-- -- --
1|1 4 7
2|2 5 8
3|3 6 9
>m.colNames();
["c1","c2","c3"]
>m.rowNames();
[1,2,3]
此外,通过 SQL 语句也可以生成矩阵:
exec
+ pivot by
: SQL 语句中,可以通过 exec 语句搭配 pivot by 子句,将表数据转换成一个矩阵形式的面板数据。更多使用场景可以参考面板数据处理教程。
>sym = `C`MS`MS`MS`IBM`IBM`C`C`C$SYMBOL
>price = 49.6 29.46 29.52 30.02 174.97 175.23 50.76 50.32 51.29
>qty = 2200 1900 2100 3200 6800 5400 1300 2500 8800 >timestamp = [09:34:07,09:35:42,09:36:51,09:36:59,09:35:47,09:36:26,09:34:16,09:35:26,09:36:12]
>t2 = table(timestamp, sym, qty, price)
>b = exec count(price) from t2 pivot by timestamp.minute(), sym;
>b;
C IBM MS
- --- --
09:34m|2
09:35m|1 1 1
09:36m|1 1 2
>typestr b;
FAST DOUBLE MATRIX
除了使用 pivot by 子句外,还可以使用以下函数生成面板数据:
pivot
:在二维维度上重组数据,其功能与 pivot by 等同。panel
:指定一个或多个指标列,生成一个或多个矩阵。
使用函数 cross
将指定函数应用于 X 和 Y 中元素的两两组合,其中 X, Y 可以是数据对向量或矩阵。
例:假设 X,Y 是向量,X 有 m 个元素,矩阵 Y 有 n 个元素,则 cross 函数将返回一个 m×n 矩阵。
>x=1 2;
>y=3 5 7;
>cross(pow, x, y);
3 5 7
- -- ---
1|1 1 1
2|8 32 128
2.2. 访问矩阵
矩阵的维度可以通过函数 shape
获得。单独获取行数和列数,可通过调用 rows
和 cols
函数实现。
>m=1..10$2:5;
>shape(m);
2 : 5
>rows(m)
2
>cols(m)
5
DolphinDB 提供了多样化的矩阵元素访问方式。
2.2.1. 根据行列索引访问矩阵
定义一个普通矩阵:
m = 1..100$10:10;
(1)获取单个单元格的值:slice
, cell
>m.slice(1,1) // 等同于 m[1,1]
12
>m.cell(1,1)
12
(2)获取多个单元格的值:slice
, cell
>m.slice([0,1],[1,2,3]) // 等同于 m[[0,1],[1,2,3]]
col1 col2 col3
11 21 31
12 22 32
>m.cells([0,1],[1,2])
[11, 22]
cells 与 slice 的区别在于:
cells 输入的行列坐标一一对应,因此返回的结果与输入的行/列坐标的长度相同;slice 获取的是行/列坐标的组合值,因此返回的结果是一个 “行坐标长度 × 列坐标长度” 维度的矩阵。
cells 返回一个向量;slice 是对矩阵进行切片,返回一个矩阵(可调用 flatten 转换为向量)。
(3)获取某列/某几列的值:slice
, col
, loc
>m.slice(0) // 等同于 m[0] 或 m[, 0]
[1,2,3,4,5,6,7,8,9,10]
>m.slice([0,3]) // 等同于 m[[0,3]]
col1 col2
1 31
2 32
3 33
4 34
5 35
6 36
7 37
8 38
9 39
10 40
>m.slice(0:3) // 等同于 m[0:3] 或 m[,0:3]
col1 col2 col3
1 11 21
2 12 22
3 13 23
4 14 24
5 15 25
6 16 26
7 17 27
8 18 28
9 19 29
10 20 30
>m.col(0)
[1,2,3,4,5,6,7,8,9,10]
>m.col(0:3)
col1 col2 col3
1 11 21
2 12 22
3 13 23
4 14 24
5 15 25
6 16 26
7 17 27
8 18 28
9 19 29
10 20 30
>m.loc(,[true, false, false, false, false, false, false, false, false, false])
col1
1
2
3
4
5
6
7
8
9
10
>m.loc(,[true, false, true, false, false, false, true, false, false, false])
col1 col2 col3
1 21 61
2 22 62
3 23 63
4 24 64
5 25 65
6 26 66
7 27 67
8 28 68
9 29 69
10 30 70
(4)获取某行/某几行的值:slice
, row
, loc
>m.slice(1,).flatten() // 等同于 m[index, ].flatten()
[2,12,22,32,42,52,62,72,82,92]
>m.slice([1,2],)
col1 col2 col3 col4 col5 col6 col7 col8 col9 col10
1 11 21 31 41 51 61 71 81 91
2 12 22 32 42 52 62 72 82 92
>m.row(1)
[2,12,22,32,42,52,62,72,82,92]
>m.row(0:2)
col1 col2 col3 col4 col5 col6 col7 col8 col9 col10
1 11 21 31 41 51 61 71 81 91
2 12 22 32 42 52 62 72 82 92
>m.loc([true, false, false, false, false, false, false, false, false, false])
col1 col2 col3 col4 col5 col6 col7 col8 col9 col10
1 11 21 31 41 51 61 71 81 91
>m.loc([true, false, true, false, false, false, false, false, false, false])
col1 col2 col3 col4 col5 col6 col7 col8 col9 col10
1 11 21 31 41 51 61 71 81 91
3 13 23 33 43 53 63 73 83 93
(5)获取指定位置的矩阵切片:slice
, loc
>m.slice(1:3,0:4) // s[1:3,0:4]
col1 col2 col3 col4
2 12 22 32
3 13 23 33
>m.slice([1,4], [2])
col1
22
25
>m.loc([true, false, true, false, false, false, false, false, false, false],[true, true, true, false, false, false, false, false, false, false])
col1 col2 col3
1 11 21
3 13 23
2.2.2. 根据标签访问矩阵
通过 rename!
函数为矩阵打上标签。
>m = 1..25$5:5;
>rlabel = 1..5
>clabel = 2022.01.01 + 1..5
>m.rename!(rlabel, clabel)
>m;
label 2022.01.02 2022.01.03 2022.01.04 2022.01.05 2022.01.06
1 1 6 11 16 21
2 2 7 12 17 22
3 3 8 13 18 23
4 4 9 14 19 24
5 5 10 15 20 25
>m.loc(1, 2022.01.03)
6
>m.loc([1,3,4], [2022.01.02, 2022.01.06])
label 2022.01.02 2022.01.06
1 1 21
3 3 23
4 4 24
// 使用数据对访问标签时,要求矩阵必须是索引矩阵或索引序列
>m.setIndexedMatrix!()
>m.loc(1:3, 2022.01.01:2022.01.03)
label 2022.01.02 2022.01.03
1 1 6
2 2 7
3 3 8
2.3. 修改矩阵
2.3.1. 追加数据
通过 append!
函数追加向量到矩阵,该向量的长度必须是矩阵行数的倍数。
>m=1..6$2:3;
>m;
#0 #1 #2
-- -- --
1 3 5
2 4 6
// 追加1列到m
>append!(m, 7 9);
#0 #1 #2 #3
-- -- -- --
1 3 5 7
2 4 6 9
// 追加2列到m
>append!(m, 8 6 1 2);
#0 #1 #2 #3 #4 #5
-- -- -- -- -- --
1 3 5 7 8 1
2 4 6 9 6 2
// 要追加的向量长度必须能被矩阵的行数除尽。
>append!(m, 3 4 5);
The size of the vector to append must be divisible by the number of rows of the matrix.
2.3.2. 修改数据
访问矩阵固定位置的元素并通过赋值修改对应元素的值。
修改列:使用 m[index] = X 修改某列;使用 m[start:end] = X 修改多列。
修改行:使用 m[index,] = X 修改某行;使用 m[start:end,] = X 修改多行。
修改矩阵窗口:使用 m[r1:r2, c1:c2] = X 来修改矩阵窗口。
其中,X 是一个标量或者向量。
2.4. 按列对矩阵进行过滤
mask
函数会过滤满足条件表达式的结果,并保留不满足条件的结果,不改变矩阵的维度。at 保留满足条件的结果,其行为和 mask 相反。
>mask(m, m%2!=0);
// 等价于 iif(m%2!=0, NULL, m)
col1 col2 col3 col4 col5 col6 col7 col8 col9 col10
2 12 22 32 42 52 62 72 82 92
4 14 24 34 44 54 64 74 84 94
6 16 26 36 46 56 66 76 86 96
8 18 28 38 48 58 68 78 88 98
10 20 30 40 50 60 70 80 90 100
>m.at(m%2!=0) // m[m%2!=0]
// 等价于 iif(m%2!=0, m, NULL)
col1 col2 col3 col4 col5 col6 col7 col8 col9 col10
1 11 21 31 41 51 61 71 81 91
3 13 23 33 43 53 63 73 83 93
5 15 25 35 45 55 65 75 85 95
7 17 27 37 47 57 67 77 87 97
9 19 29 39 49 59 69 79 89 99
使用 lambda 表达式对矩阵的每一个列进行过滤。注意,按列对矩阵进行过滤时,lambda 表达式只能接受一个参数,并且返回的结果必须是 BOOL 类型的标量。
>m=matrix(0 2 3 4,0 0 0 0,4 7 8 2);
//返回矩阵中不全为0的列
>m[x->!all(x==0)];
#0 #1
-- --
0 4
2 7
3 8
4 2
//返回矩阵中均值大于4的列
>m=matrix(0 2 3 4,5 3 6 9,4 7 8 2);
>m[def (x):avg(x)>4];
#0 #1
-- --
5 4
3 7
6 8
9 2
通过 iif
函数修改矩阵中满足条件的元素。
>m1=1..6$3:2
>m2=6..1$3:2
>iif(m1>m2, m1, m2);
col1 col2
6 4
5 5
4 6
>m=matrix(NULL 2 3 4,NULL NULL 3 NULL,4 7 8 2);
>iif(isNull(m), 0, m)
col1 col2 col3
0 0 4
2 0 7
3 3 8
4 0 2
2.5. 矩阵的拼接
通过函数 concatMatrix
将多个矩阵进行水平/垂直的拼接:
>m1 = matrix(4 0 5, 2 1 8);
>m2 = matrix(2 9 8, 3 7 -3, 6 4 2, 0 5 8);
>m3 = matrix(1 -1 6 2, 1 -3 1 9, 5 3 0 -4, 1 NULL 3 4);
>concatMatrix([m1, m2]);
col1 col2 col3 col4 col5 col6
4 2 2 3 6 0
0 1 9 7 4 5
5 8 8 (3) 2 8
>concatMatrix([m2, m3], false);
col1 col2 col3 col4
2 3 6 0
9 7 4 5
8 (3) 2 8
1 1 5 1
(1) (3) 3
6 1 0 3
2 9 (4) 4
2.6. 矩阵的连接和对齐
矩阵的连接和表类似,可以根据行标签进行 'inner join', 'outer join', 'left join', 'right join', 与'asof join' 等连接操作,合并两个矩阵。通过函数 merge
实现。
>m1 = matrix([1.2, 7.8, 4.6, 5.1, 9.5], [0.15, 1.26, 0.45, 1.02, 0.33]).rename!([2012.01.01, 2015.02.01, 2015.03.01, 2015.04.01, 2015.05.01], `x1`x2).setIndexedMatrix!()
>m2 = matrix([1.0, 2.0, 3.0, 4.0], [0.14, 0.26, 0.35, 0.48]).rename!([2015.02.01, 2015.02.16, 2015.05.01, 2015.05.02], `y1`y2).setIndexedMatrix!()
>m1;
label x1 x2 y1 y2
2012.01.01 1.2 0.15
2015.02.01 7.8 1.26 1 0.14
2015.03.01 4.6 0.45 2 0.26
2015.04.01 5.1 1.02 2 0.26
2015.05.01 9.5 0.33 3 0.35
>m2;
label y1 y2
2015.02.01 1 0.14
2015.02.16 2 0.26
2015.05.01 3 0.35
2015.05.02 4 0.48
>merge(m1, m2, 'asof');
label x1 x2 y1 y2
2012.01.01 1.2 0.15
2015.02.01 7.8 1.26 1 0.14
2015.03.01 4.6 0.45 2 0.26
2015.04.01 5.1 1.02 2 0.26
2015.05.01 9.5 0.33 3 0.35
>merge(m1, m2, 'outer');
label x1 x2 y1 y2
2012.01.01 1.2 0.15
2015.02.01 7.8 1.26 1 0.14
2015.02.16 2 0.26
2015.03.01 4.6 0.45
2015.04.01 5.1 1.02
2015.05.01 9.5 0.33 3 0.35
2015.05.02 4 0.48
连接操作除了用于合并矩阵外,还可用于矩阵标签和数据的对齐。通过函数 align
实现。
align 支持按行对齐,按列对齐或者同时按行列进行矩阵对齐,详情请参考用户手册的说明。
>x1 = [09:00:00, 09:00:01, 09:00:03]
>x2 = [09:00:00, 09:00:02, 09:00:03, 09:00:04]
>m1 = matrix(1 2 3, 2 3 4, 3 4 5).rename!(x1)
>m2 = matrix(1 2 3, 2 3 4, 3 4 5, 4 5 6).rename!(x2)
>m = align(m1, m3, 'fj', false);
>m[0];
09:00:00 09:00:01 09:00:03 09:00:04
1 2 3
2 3 4
3 4 5
>m[1];
09:00:00 09:00:01 09:00:03 09:00:04
1 2 3
2 3 4
3 4 5
4 5 6
2.7. 矩阵空值处理
2.7.1. 填充空值
矩阵的二元计算(非聚合运算)中,若包含空值,返回的结果也为 NULL。参考手册第 6 章:NULL 值的操作。
>m1 = 3 -2 1 NULL -5 6$2:3
col1 col2 col3
3 1 (5)
(2) 6
>m2 = 5 NULL 2 4 -5 9$2:3
col1 col2 col3
5 2 (5)
4 9
>m1 + m2
col1 col2 col3
8 3 (10)
15
可以看到,包含空值的单元格计算结果也为空,这可能与预期的结果不符。某些场景下,用户希望空值不影响计算结果,或者以某个特定的值去填充空值。
常规的空值填充方法有:
bfill
和bfill!
:使用 NULL 后的非NULL元素填充 NULL 值。ffill
和ffill!
:使用 NULL 前的非NULL元素填充 NULL 值。lfill
和lfill!
:线性填充非 NULL 元素之间的 NULL 值。fill!
,nullFill
和nullFill!
:用指定值填充 NULL 值。
上述方法都需要对参与计算的矩阵单独填充。对于矩阵间的二元计算,可以通过 withNullFill
函数同时实现填充和计算。
2.7.2. 去除空行/列
若想删除矩阵中全为空值或空值较多的行/列,可以通过 dropna
实现。
下例使用 dropna 删除 mask 过滤后全为空的行。
>m=1..100$10:10
>dropna(m.at(m%2!=0))
col1 col2 col3 col4 col5 col6 col7 col8 col9 col10
1 11 21 31 41 51 61 71 81 91
3 13 23 33 43 53 63 73 83 93
5 15 25 35 45 55 65 75 85 95
7 17 27 37 47 57 67 77 87 97
9 19 29 39 49 59 69 79 89 99
3. 矩阵的运算
3.1. 矩阵基本运算
DolphinDB 中矩阵的每一列都可以视为一个向量,所以大部分向量函数和聚合函数都适用于矩阵,计算结果等价于对矩阵每列进行计算后结果的拼接。例如:
>m=1..6$2:3;
>m;
#0 #1 #2
-- -- --
1 3 5
2 4 6
// 每个元素的余弦值
>cos(m);
#0 #1 #2
--------- --------- --------
0.540302 -0.989992 0.283662
-0.416147 -0.653644 0.96017
DolphinDB 中矩阵可以直接与其它矩阵或标量进行四则运算。
矩阵和标量进行运算,等价于矩阵每个元素和标量进行运算。注意若标量为 NULL,计算结果也为 NULL:
>m=1..10$5:2;
>m;
#0 #1
-- --
1 6
2 7
3 8
4 9
5 10
>2.1*m;
#0 #1
---- ----
2.1 12.6
4.2 14.7
6.3 16.8
8.4 18.9
10.5 21
>m\2;
#0 #1
--- ---
0.5 3
1 3.5
1.5 4
2 4.5
2.5 5
>m+1.1;
#0 #1
--- ----
2.1 7.1
3.1 8.1
4.1 9.1
5.1 10.1
6.1 11.1
>m*NULL;
#0 #1
-- --
矩阵和矩阵进行运算时,参与计算的两个矩阵的维度必须一致:
>m=1..10$5:2;
>n=3..12$5:2;
>m;
#0 #1
-- --
1 6
2 7
3 8
4 9
5 10
>n;
#0 #1
-- --
3 8
4 9
5 10
6 11
7 12
>m*n;
#0 #1
-- ---
3 48
8 63
15 80
24 99
35 120
>m\n;
#0 #1
-------- --------
0.333333 0.75
0.5 0.777778
0.6 0.8
0.666667 0.818182
0.714286 0.833333
>m+n;
#0 #1
-- --
4 14
6 16
8 18
10 20
12 22
>n=3..10$4:2
>m+n;
Incompatible vector size
DolphinDB 中部分矩阵运算,如相乘、求逆、求行列式、矩阵分解等采用了 OpenBlas 和 Lapack 进行优化,与 Matlab 的性能在一个数量级。
3.2. 矩阵相乘,求逆,转置,求矩阵行列式
矩阵相乘:
dot(X,Y) 或 X**Y:返回 X 和 Y 的矩阵乘法的结果。
>x=1..6$2:3;
>y=1 2 3;
>x;
#0 #1 #2
-- -- --
1 3 5
2 4 6
>y;
[1,2,3]
>x dot y;
#0
--
22
28
>y=6..1$3:2;
>y;
#0 #1
-- --
6 3
5 2
4 1
>x**y;
#0 #1
-- --
41 14
56 20
>y**x;
#0 #1 #2
-- -- --
12 30 48
9 23 37
6 16 26
矩阵求逆:
inverse(X):如果 X 可逆,返回矩阵 X 的逆矩阵。
>x=1..4$2:2;
>x;
#0 #1
-- --
1 3
2 4
>x.inverse();
#0 #1
-- --
-2 1.5
1 -0.5
矩阵转置:
transpose(X):如果 X 为矩阵,返回 X 的转置矩阵。
>x=1..6 $ 3:2;
>x;
#0 #1
-- --
1 4
2 5
3 6
>transpose x;
#0 #1 #2
-- -- --
1 2 3
4 5 6
求行列式:
det(X):计算矩阵的行列式,在计算中,NULL 值用 0 代替。
>x=1..4$2:2;
>x;
#0 #1
-- --
1 3
2 4
>X.det();
-2
>x=1 2 3 6 5 4 8 7 NULL $3:3;
>x;
#0 #1 #2
-- -- --
1 6 8
2 5 7
3 4
>det x;
42
>x=1 2 3 6 5 4 3 6 9$3:3;
>x;
#0 #1 #2
-- -- --
1 6 3
2 5 6
3 4 9
>det x;
0
3.3. 按行计算和按列计算
“矩阵基本运算”一节简单介绍了函数在矩阵上应用时的计算规则,可以看出对矩阵直接应用向量函数和聚合函数,都是按列进行计算的。此外,DolphinDB 手册提供了部分按行计算的函数和对应的高阶函数byRow
。
部分函数不支持矩阵的列计算,此时可以采用高阶函数 each
实现:
//两个矩阵按列进行矩阵相乘运算
>x=1..6$2:3;
>y=6..1$2:3;
>x;
#0 #1 #2
-- -- --
1 3 5
2 4 6
>y;
#0 #1 #2
-- -- --
6 4 2
5 3 1
>each(**, x, y);
[16,24,16] //比如,24=3*4+4*3
4. 矩阵分解与求解线性方程
DolphinDB 实现了以下矩阵分解函数:
4.1. lu 分解
lu
:对矩阵进行 lu分 解,X = P**L**U,其中P为置换矩阵,L 为下三角矩阵,U 为上三角矩阵。如果 permute=true,返回 PL=P**L 和 U。
>m=matrix([2 5 7 5, 5 2 5 4, 8 2 6 4, 7 8 6 8]);
>m;
#0 #1 #2 #3
-- -- -- --
2 5 8 7
5 2 2 8
7 5 6 6
5 4 4 8
>p,l,u=lu(m);
>p;
#0 #1 #2 #3
-- -- -- --
0 1 0 0
0 0 0 1
1 0 0 0
0 0 1 0
>l;
#0 #1 #2 #3
-------- ----- --------- --
1 0 0 0
0.285714 1 0 0
0.714286 0.12 1 0
0.714286 -0.44 -0.461538 1
>u;
#0 #1 #2 #3
-- -------- -------- --------
7 5 6 6
0 3.571429 6.285714 5.285714
0 0 -1.04 3.08
0 0 0 7.461538
>pl,u=lu(m,true);
>pl;
#0 #1 #2 #3
-------- ----- --------- --
0.285714 1 0 0
0.714286 -0.44 -0.461538 1
1 0 0 0
0.714286 0.12 1 0
>u;
#0 #1 #2 #3
-- -------- -------- --------
7 5 6 6
0 3.571429 6.285714 5.285714
0 0 -1.04 3.08
0 0 0 7.461538
可以对非方阵的矩阵进行 lu 分解。
>m=matrix([2 5 5, 5 5 4, 8 6 4, 7 6 8]);
>m;
#0 #1 #2 #3
-- -- -- --
2 5 8 7
5 5 6 6
5 4 4 8
>p,l,u=lu(m);
>p;
#0 #1 #2
-- -- --
0 1 0
1 0 0
0 0 1
>l;
#0 #1 #2
--- --------- --
1 0 0
0.4 1 0
1 -0.333333 1
>u;
#0 #1 #2 #3
-- -- --------- --------
5 5 6 6
0 3 5.6 4.6
0 0 -0.133333 3.533333
>pl,u=lu(m,true);
>pl;
#0 #1 #2
--- --------- --
0.4 1 0
1 0 0
1 -0.333333 1
>u;
#0 #1 #2 #3
-- -- --------- --------
5 5 6 6
0 3 5.6 4.6
0 0 -0.133333 3.533333
4.2. qr 分解
qr
:X 为一个 m*n 的矩阵,对 X 进行 qr 分解,X=Q**R,Q 为正交矩阵,R 为上三角矩阵。mode 的选项有:'full', 'r', 'economic', 'raw'。
mode 为'full':返回完整的 qr 分解结果,即 Q 为 m*m 的矩阵,R 为 m*n 的矩阵。
>m=matrix([2 5 5, 5 5 4, 8 6 4, 7 6 8]);
>m;
#0 #1 #2 #3
-- -- -- --
2 5 8 7
5 5 6 6
5 4 4 8
>q,r=qr(m,mode='full',pivoting=false);
>q;
#0 #1 #2
--------- --------- ---------
-0.272166 0.93784 0.215365
-0.680414 -0.029307 -0.732242
-0.680414 -0.345828 0.646096
>r;
#0 #1 #2 #3
--------- --------- --------- ----------
-7.348469 -7.484552 -8.981462 -11.430952
0 3.159348 5.943561 3.622407
0 0 -0.086146 2.282872
>m=matrix([2 5 7 5, 5 2 5 4, 8 2 6 4]);
>m;
#0 #1 #2
-- -- --
2 5 8
5 2 2
7 5 6
5 4 4
>q,r=qr(m); //mode='full',pivoting=false
>q;
#0 #1 #2 #3
--------- --------- --------- ---------
-0.197066 0.903357 0.300275 0.234404
-0.492665 -0.418267 0.459245 0.609449
-0.68973 -0.02475 0.170745 -0.703211
-0.492665 0.091573 -0.818398 0.281284
>r;
#0 #1 #2
---------- -------- ---------
-10.148892 -7.38997 -8.670898
0 3.922799 6.608121
0 0 1.071571
0 0 0
对于 m>n 的矩阵,mode 为 full 时 R 矩阵的下面的 m−n 行全为 0,X=Q**R 可以进一步分解为: X = Q**([R1,0]T) = [Q1,Q2]**([R1,0]T) = Q1**R1。
mode 取 'economic': m>n,即返回 Q1 和 R1,Q1 为 m*n 的矩阵,R 为 n*n 的矩阵;m<=n 时,返回结果和 mode 取 'full' 时一致。
>m=matrix([2 5 5, 5 5 4, 8 6 4, 7 6 8]);
>m;
#0 #1 #2 #3
-- -- -- --
2 5 8 7
5 5 6 6
5 4 4 8
>q,r=qr(m,mode='economic',pivoting=false);
>q;
#0 #1 #2
--------- --------- ---------
-0.272166 0.93784 0.215365
-0.680414 -0.029307 -0.732242
-0.680414 -0.345828 0.646096
>r;
#0 #1 #2 #3
--------- --------- --------- ----------
-7.348469 -7.484552 -8.981462 -11.430952
0 3.159348 5.943561 3.622407
0 0 -0.086146 2.282872
>m=matrix([2 5 7 5, 5 2 5 4, 8 2 6 4]);
>m;
#0 #1 #2
-- -- --
2 5 8
5 2 2
7 5 6
5 4 4
>q,r=qr(m,mode='economic'); //pivoting=false
>q;
#0 #1 #2
--------- --------- ---------
-0.197066 0.903357 0.300275
-0.492665 -0.418267 0.459245
-0.68973 -0.02475 0.170745
-0.492665 0.091573 -0.818398
>r;
#0 #1 #2
---------- -------- ---------
-10.148892 -7.38997 -8.670898
0 3.922799 6.608121
0 0 1.071571
mode取 'r':返回 qr(X,mode = 'full') 中的R。
>m=matrix([2 5 5, 5 5 4, 8 6 4, 7 6 8]);
>m;
#0 #1 #2 #3
-- -- -- --
2 5 8 7
5 5 6 6
5 4 4 8
>r=qr(m,mode='r',pivoting=false);
>r;
#0 #1 #2 #3
--------- --------- --------- ----------
-7.348469 -7.484552 -8.981462 -11.430952
0 3.159348 5.943561 3.622407
0 0 -0.086146 2.282872
>m=matrix([2 5 7 5, 5 2 5 4, 8 2 6 4]);
>m;
#0 #1 #2
-- -- --
2 5 8
5 2 2
7 5 6
5 4 4
>r=qr(m,mode='r'); //pivoting=false
>r;
#0 #1 #2
---------- -------- ---------
-10.148892 -7.38997 -8.670898
0 3.922799 6.608121
0 0 1.071571
0 0 0
mode 取 'raw': 返回矩阵 h、数组 tau 和矩阵 R。qr 分解的计算通过 householder 变换实现,h 为其计算 R 矩阵时的变换矩阵,tau 为 householder 变换的变换因子。
>m=matrix([2 5 5, 5 5 4, 8 6 4, 7 6 8]);
>m;
#0 #1 #2 #3
-- -- -- --
2 5 8 7
5 5 6 6
5 4 4 8
>h,tau,r=qr(m,mode=`raw,pivoting=false);
>h;
#0 #1 #2 #3
--------- --------- --------- ----------
-7.348469 -7.484552 -8.981462 -11.430952
0.534847 3.159348 5.943561 3.622407
0.534847 0.553547 -0.086146 2.282872
>tau;
[1.272166,1.530908,0]
>r;
#0 #1 #2 #3
--------- --------- --------- ----------
-7.348469 -7.484552 -8.981462 -11.430952
0 3.159348 5.943561 3.622407
0 0 -0.086146 2.282872
pivoting 取 true:计算秩显 qr 分解,计算 A**P=Q**R,P 是置换矩阵,R 矩阵满足对角线上的元素不递减。qr(X,pivoting=true,[mode='full']) 返回矩阵 Q,R 和 piv 数组,piv 为 P 矩阵的置换规则。
>m=matrix([2 5 5, 5 5 4, 8 6 4, 7 6 8]);
>m;
#0 #1 #2 #3
-- -- -- --
2 5 8 7
5 5 6 6
5 4 4 8
>q,r,piv=qr(m,mode=`full,pivoting=true);
>q;
#0 #1 #2
--------- --------- --------
-0.742781 0.629178 0.228934
-0.557086 -0.391111 -0.73259
-0.371391 -0.67169 0.641016
>r;
#0 #1 #2 #3
--------- --------- ---------- ---------
-10.77033 -6.127946 -11.513111 -7.9849
0 -4.055647 -3.315938 -1.496423
0 0 2.33513 0.045787
>piv;
[2,0,3,1]
//置换规则: 如果 piv[i] != i, 则m[i] = m[piv[i]]
>q**r;
#0 #1 #2 #3
-- -- -- --
8 2 7 5
6 5 6 5
4 5 8 4
4.3. svd 分解
函数 svd
计算矩阵 svd 分解:X = U ** Σ ** VT,X 为 m*n 的矩阵,输出 U, s(diag(Σ)), Vt(V.transpose());U 和 V 均为单位正交阵,U 为左奇异矩阵,V 为右奇异矩阵。Σ 仅在主对角线上有值,s 为矩阵的奇异值。
s 是长度为 k 的一维矩阵;fullfullMatrice 为 true 时,U 和 Vh 的维度分别为(m, m) 与 (n, n);fullfullMatrice 为 false 时,U 和 Vh 的维度分别为(m, k), (k, n), k = min(m, n)。 computeUV 为 false 时,只返回 s。
>m=matrix([2 5 5, 5 5 4, 8 6 4, 7 6 8]);
>m;
#0 #1 #2 #3
-- -- -- --
2 5 8 7
5 5 6 6
5 4 4 8
>u,s,vh = svd(m); //fullMatrices=true,computeUV=true
>u;
#0 #1 #2
--------- --------- ---------
-0.604057 -0.734356 0.309574
-0.570062 0.126705 -0.811773
-0.556906 0.666834 0.495165
>s;
[19.202978,3.644527,1.721349]
>vh;
#0 #1 #2 #3
--------- --------- --------- ---------
-0.356349 -0.421717 -0.545772 -0.63032
0.68568 -0.101776 -0.671497 0.261873
-0.559964 -0.308094 -0.240155 0.730646
-0.29883 0.846684 -0.439944 -0.016602
>u,s,vh = svd(m, fullMatrices=false); //computeUV=true
>u;
#0 #1 #2
--------- --------- ---------
-0.604057 -0.734356 0.309574
-0.570062 0.126705 -0.811773
-0.556906 0.666834 0.495165
>s;
[19.202978,3.644527,1.721349]
>vh;
#0 #1 #2 #3
--------- --------- --------- --------
-0.356349 -0.421717 -0.545772 -0.63032
0.68568 -0.101776 -0.671497 0.261873
-0.559964 -0.308094 -0.240155 0.730646
>s = svd(m, computeUV=false);
>s;
[19.202978,3.644527,1.721349]
4.4. cholesky 分解
cholesky(X, [lower=true]): 对矩阵进行 Cholesky 分解 X = L ** LT, X 是一个对称正定矩阵; lower 是一个布尔值,表示是否使用输入矩阵的下三角来计算分解。默认值为 true,表示使用下三角计算。如果 lower 为 false,表示使用上三角计算。
>m=[1, 0, 1, 0, 2, 0, 1, 0, 3]$3:3;
>m;
#0 #1 #2
-- -- --
1 0 1
0 2 0
1 0 3
>L=cholesky(m);
>L;
#0 #1 #2
-- -------- --------
1 0 0
0 1.414214 0
1 0 1.414214
>L**transpose(L);
#0 #1 #2
-- -- --
1 0 1
0 2 0
1 0 3
>cholesky(m, false);
#0 #1 #2
-- -------- --------
1 0 1
0 1.414214 0
0 0 1.414214
4.5. schur 分解
schur
: 对矩阵进行 schur 分解,X = U ** T ** UH;U 为幺正矩阵 (U-1 = UT),T 为上三角矩阵,T 的对角元素是 A 的所有特征值;sort 根据指定的特征值顺序对分解得到的矩阵进行重新排序,默认为 NULL,目前支持的有 'lhp': 左半平面 (e < 0.0);'rhp': 右半平面 (e > 0.0);'iuc': 单位圆盘的内部 (abs(e) < 1);'ouc': 单位圆盘的外部 (abs(e) >= 1)。
>m=matrix([2 5 7 5, 5 2 5 4, 8 2 6 4, 7 8 6 8]);
>m;
#0 #1 #2 #3
-- -- -- --
2 5 8 7
5 2 2 8
7 5 6 6
5 4 4 8
>t,u=schur(m);
>t;
#0 #1 #2 #3
-------- --------- --------- ---------
21.16354 -1.073588 -0.473548 -4.270044
0 -4.306007 -1.391659 2.039609
0 0 -0.995651 -2.879786
0 0 0 2.138117
>u;
#0 #1 #2 #3
-------- --------- --------- ---------
0.52214 0.818236 0.198151 0.136364
0.401387 -0.461653 0.785756 0.091394
0.568479 -0.320408 -0.540423 0.531143
0.493041 -0.121263 -0.226421 -0.831228
>u**t**u.transpose();
#0 #1 #2 #3
-- -- -- --
2 5 8 7
5 2 2 8
7 5 6 6
5 4 4 8
>t,u, sdim=schur(m,'lhp') //'lhp':(e < 0.0), sdim 为符合sort条件的特征值数量
>t;
#0 #1 #2 #3
--------- --------- --------- ---------
-4.306007 -1.390041 -1.099778 -1.857969
0 -0.995651 -0.414519 2.960681
0 0 21.16354 -4.29753
0 0 0 2.138117
>u;
#0 #1 #2 #3
-------- --------- -------- ---------
-0.8395 -0.207229 0.483426 0.136364
0.444339 -0.793483 0.405703 0.091394
0.296182 0.529453 0.591475 0.531143
0.100391 0.217073 0.501858 -0.831228
>sdim;
2
>t,u, sdim=schur(m,'rhp') //'rhp':(e < 0.0)
>t;
#0 #1 #2 #3
-------- --------- --------- ---------
21.16354 -3.020883 0.002732 3.237958
0 2.138117 2.443445 -2.818711
0 0 -4.306007 -0.688714
0 0 0 -0.995651
>u;
#0 #1 #2 #3
-------- --------- --------- ---------
0.52214 0.258617 -0.777021 -0.238171
0.401387 -0.597888 0.267022 -0.640405
0.568479 0.594002 0.567898 0.03853
0.493041 -0.472026 -0.049293 0.729158
>sdim;
2
>t,u, sdim=schur(m,'iuc') //'iuc':(abs(e) < 1.0)
>t;
#0 #1 #2 #3
--------- -------- --------- ---------
-0.995651 -0.02048 1.390574 3.449116
0 21.16354 1.174494 -4.266858
0 0 -4.306007 -0.764178
0 0 0 2.138117
>u;
#0 #1 #2 #3
--------- -------- --------- ---------
0.133953 0.522264 -0.831085 0.136364
-0.903631 0.400552 0.121062 0.091394
0.373493 0.568825 0.504805 0.531143
0.161277 0.49319 0.199534 -0.831228
>sdim;
1
>t,u, sdim=schur(m,'ouc') //'ouc':(abs(e) >= 1.0)
>t;
#0 #1 #2 #3
-------- --------- --------- ---------
21.16354 -1.073588 -2.823677 3.237958
0 -4.306007 2.443445 -0.355379
0 0 2.138117 -2.879786
0 0 0 -0.995651
>u;
#0 #1 #2 #3
-------- --------- --------- ---------
0.52214 0.818236 -0.03367 -0.238171
0.401387 -0.461653 -0.464378 -0.640405
0.568479 -0.320408 0.75676 0.03853
0.493041 -0.121263 -0.45884 0.729158
0.161277 0.49319 0.199534 -0.831228
>sdim;
3
4.6. 求解线性方程组
solve(a,y): 求解线性方程组 a*x=y
>a=[1,0,2,1,2,5,1,5,-1]$3:3;
>y=[6,-4,27];
>a;
#0 #1 #2
-- -- --
1 1 1
0 2 5
2 5 -1
>x = solve(a,y);
>x;
[5,3,-2]
>a ** matrix(x);
#0
--
6
-4
27
5. 矩阵高级运算
5.1. 计算矩阵特征值和特征向量
函数 eig(X)
计算矩阵的特征值和特征向量。返回一个字典,包含两个键:values(特征值)与 vectors(特征向量)。
>m=matrix([2 5 7 5, 5 2 5 4, 8 2 6 4, 7 8 6 8]);
>m;
#0 #1 #2 #3
-- -- -- --
2 5 8 7
5 2 2 8
7 5 6 6
5 4 4 8
>v=eig(m);
>v[`values];
[19.750181,-3.807927,-1.551136,3.608882]
>v[`vectors];
#0 #1 #2 #3
-------- --------- --------- ---------
0.485606 -0.853982 -0.034777 -0.183556
0.413406 0.301775 -0.845881 -0.15004
0.553396 0.40595 0.507665 -0.520802
0.535757 0.121868 0.15985 0.820098
5.2. PCA 计算
pca
: 对数据源中指定列中的数据求主成分分析。
返回的结果是一个字典,包含以下键:
- explainedVarianceRatio: 对应一个长度为k的向量,包含前k个主成分分别解释的方差权重。
- singularValues: 对应一个长度为k的向量,包含主成分方差(协方差矩阵特征值)。
- components: 对应长度为size(xColNames)*k的主成分分析矩阵。
>x = [7,1,1,0,5,2]
>y = [0.7, 0.9, 0.01, 0.8, 0.09, 0.23]
>t=table(x, y)
>ds = sqlDS(<select * from t>);
>pca(ds);
components->
#0 #1
--------- ---------
-0.999883 0.015306
-0.015306 -0.999883
explainedVarianceRatio->[0.980301,0.019699]
singularValues->[6.110802,0.866243]
6. JIT 的支持
从 1.20 版本开始,DolphinDB database 的 JIT 增加了对矩阵运算的支持。支持矩阵作为函数参数和返回值,支持矩阵的四则运算,函数 det
与 flatten
,以及矩阵的转置、点乘等运算。
//定义对角矩阵,jit计算比非jit快了10倍左右
@jit
def diagonalMatrix_jit(data, m){
n=data.size()
res=m
for( i in 0:n){
//i in 0:n
res[i*n+i]=data[i]
}
return res
}
def diagonalMatrix(data, m){
n=data.size()
res=m
for( i in 0:n){
x=m[i]
for(j in 0:n){
if(i==j){
x[j]=data[i]
}
}
res[i]=x
}
return res
}
>m = matrix(DOUBLE,32,32)
>data=1..32
>timer(1000) diagonalMatrix(data,m)
Time elapsed: 420.021 ms
>timer(1000) diagonalMatrix_jit(data,m)
Time elapsed: 41.026 ms
//获取矩阵的上三角矩阵
@jit
def upperTriangularMatrix_jit(m, rowNum,colNum){
upperM=m
for( i in 0:colNum){
for(j in 0:rowNum){
if(i<j){
upperM[i*rowNum+j]=0
}
}
}
return upperM
}
def upperTriangularMatrix(m, rowNum,colNum){
upperM=m
for( i in 0:colNum){
col=flatten(m[i])
for(j in 0:rowNum){
if(i<j){
col[j]=0
}
}
upperM[i]=col
}
return upperM
}
>m=1..1024$32:32
timer(1000) upperTriangularMatrix_jit(m,32,32)
Time elapsed: 24.049 ms
>timer(1000) upperTriangularMatrix(m,32,32)
Time elapsed: 355.052 ms
对于矩阵的转置、点乘等函数,由于函数实现已做优化,所以 JIT 和非 JIT 的性能差别不大,其实现的目的是方便用户在 JIT 函数中使用。