DolphinDB 教程：矩阵运算

在处理面板数据时，矩阵运算相对于表运算而言，速度更快，效率更高。DolphinDB 针对矩阵运算场景为用户提供了丰富且便捷的内置函数，其中部分矩阵运算采用 OpenBlas 和 Lapack 进行优化，与 Matlab 的性能在一个数量级。本教程将重点介绍 DolphinDB 矩阵的特点及矩阵运算的相关方法。

本教程将讲述以下有关矩阵内容：

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 ** Σ ** V^T，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 ** L^T, 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 ** U^H；U 为幺正矩阵 (U^-1 = U^T)，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 函数中使用。

浙江智臾科技有限公司 / Tutorials_CN

DolphinDB 教程：矩阵运算

1. 矩阵的存储

2. 矩阵的基础操作

2.1. 矩阵的创建

2.2. 访问矩阵

2.2.1. 根据行列索引访问矩阵

2.2.2. 根据标签访问矩阵

2.3. 修改矩阵

2.3.1. 追加数据

2.3.2. 修改数据

2.4. 按列对矩阵进行过滤

2.5. 矩阵的拼接

2.6. 矩阵的连接和对齐

2.7. 矩阵空值处理

2.7.1. 填充空值

2.7.2. 去除空行/列

3. 矩阵的运算

3.1. 矩阵基本运算

3.2. 矩阵相乘，求逆，转置，求矩阵行列式

3.3. 按行计算和按列计算

4. 矩阵分解与求解线性方程

4.1. lu 分解

4.2. qr 分解

4.3. svd 分解

4.4. cholesky 分解

4.5. schur 分解

4.6. 求解线性方程组

5. 矩阵高级运算

5.1. 计算矩阵特征值和特征向量

5.2. PCA 计算

6. JIT 的支持

About

Releases

Contributors

Activities

浙江智臾科技有限公司 / Tutorials_CN .gitee-modal { width: 500px !important; }

DolphinDB 教程：矩阵运算

1. 矩阵的存储

2. 矩阵的基础操作

2.1. 矩阵的创建

2.2. 访问矩阵

2.2.1. 根据行列索引访问矩阵

2.2.2. 根据标签访问矩阵

2.3. 修改矩阵

2.3.1. 追加数据

2.3.2. 修改数据

2.4. 按列对矩阵进行过滤

2.5. 矩阵的拼接

2.6. 矩阵的连接和对齐

2.7. 矩阵空值处理

2.7.1. 填充空值

2.7.2. 去除空行/列

3. 矩阵的运算

3.1. 矩阵基本运算

3.2. 矩阵相乘，求逆，转置，求矩阵行列式

3.3. 按行计算和按列计算

4. 矩阵分解与求解线性方程

4.1. lu 分解

4.2. qr 分解

4.3. svd 分解

4.4. cholesky 分解

4.5. schur 分解

4.6. 求解线性方程组

5. 矩阵高级运算

5.1. 计算矩阵特征值和特征向量

5.2. PCA 计算

6. JIT 的支持

About

Releases

The Open Source Evaluation Index is derived from the OSS Compass evaluation system, which evaluates projects around the following three dimensions

Contributors

Activities

Search

浙江智臾科技有限公司 / Tutorials_CN