矩阵运算

在处理面板数据时，矩阵运算相对于表运算而言，速度更快，效率更高。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 ** Σ ** 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 函数中使用。