矩阵运算

在处理面板数据时,矩阵运算相对于表运算而言,速度更快,效率更高。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 获得。单独获取行数和列数,可通过调用 rowscols 函数实现。

>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

可以看到,包含空值的单元格计算结果也为空,这可能与预期的结果不符。某些场景下,用户希望空值不影响计算结果,或者以某个特定的值去填充空值。

常规的空值填充方法有:

  • bfillbfill!:使用 NULL 后的非NULL元素填充 NULL 值。

  • ffillffill!:使用 NULL 前的非NULL元素填充 NULL 值。

  • lfilllfill!:线性填充非 NULL 元素之间的 NULL 值。

  • fill!nullFillnullFill!:用指定值填充 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 增加了对矩阵运算的支持。支持矩阵作为函数参数和返回值,支持矩阵的四则运算,函数 detflatten,以及矩阵的转置、点乘等运算。

//定义对角矩阵,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 函数中使用。