DolphinDB实现了一系列矩阵运算功能，例如矩阵的创建和访问，矩阵相乘、求逆、转置，矩阵的分解，计算矩阵的特征值等。本篇教程介绍在DolphinDB中如何进行矩阵运算。本文的所有例子都基于DolphinDB 1.20。
本教程将讲述以下有关矩阵内容：
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的矩阵为column major，读取或构建矩阵时，按列操作比按行操作要更为高效得多。
给定一个长度为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毫秒。
用函数matrix创建一个指定数据类型和维度的矩阵。
// 创建一个整数矩阵，所有值都设为默认值0。
>matrix(int, 2, 3);
#0 #1 #2
-- -- --
0 0 0
0 0 0
函数matrix也可根据已有的数据包括向量、矩阵、表、元组及这些数据结构的组合创建一个矩阵。
>matrix(1 2 3);//矩阵按列存储，创建时也按列创建。
#0
--
1
2
3
>matrix([1],[2])
#0 #1
-- --
1 2
>matrix(1 2 3, 4 5 6);
#0 #1
-- --
1 4
2 5
3 6
>matrix(table(1 2 3 as id, 4 5 6 as value)); //根据table创建矩阵
#0 #1
-- --
1 4
2 5
3 6
//根据向量、矩阵与表的组合创建一个矩阵
>matrix([1 2 3, 4 5 6], 7 8 9, table(0.5 0.6 0.7 as id), 1..9$3:3);
#0 #1 #2 #3 #4 #5 #6
-- -- -- --- -- -- --
1 4 7 0.5 1 4 7
2 5 8 0.6 2 5 8
3 6 9 0.7 3 6 9
SQL中exec语句和pivot by一起使用时将table转换成一个矩阵。
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(func, funcArgs, rowAlignCol, colAlignCol)函数在二维维度上重组数据，将其转换为一个新的矩阵。
syms=`BIDU`MSFT`ORCL$SYMBOL
sym=syms[0 0 0 0 0 0 0 1 1 1 1 1 1 1 2 2 2 2 2 2 2]
price=172.12 170.32 172.25 172.55 175.1 174.85 174.5 36.45 36.15 36.3 35.9 36.5 37.15 36.9 40.1 40.2 40.25 40.15 40.1 40.05 39.95
qty=100* 10 3 7 8 25 6 10 4 5 1 2 8 6 10 2 2 5 5 4 4 3
trade_time=09:40:00+1 30 65 90 130 185 195 10 40 90 140 160 190 200 5 45 80 140 170 190 210
t1=table(sym, price, qty, trade_time);
stockprice=pivot(wavg, [t1.price, t1.qty], minute(t1.trade_time), t1.sym);
>stockprice;
BIDU MSFT ORCL
---------- --------- ---------
09:40m|171.704615 36.283333 40.15
09:41m|172.41 36.3 40.25
09:42m|175.1 36.38 40.127778
09:43m|174.63125 36.99375 40.007143
>typestr(stockprice);
FAST DOUBLE MATRIX
使用函数cross(func, X, [Y])将指定函数应用于X和Y中元素的两两组合。假设矩阵X有m列，矩阵Y有n列，如果func(X[i], Y[j])是标量，将返回一个m×n矩阵。
>x=1 2;
>y=3 5 7;
>cross(pow, x, y);
3 5 7
- -- ---
1|1 1 1
2|8 32 128
运算符cast($)可以把一个向量转换为一个矩阵。
>m=1..10$5:2;
>m;
#0 #1
-- --
1 6
2 7
3 8
4 9
5 10
>cast(m,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]
函数reshape可进行矩阵与向量之间的转换。
>m1=1..6$3:2;
>m2 = m1$2:3;
>m2;
#0 #1 #2
-- -- --
1 3 5
2 4 6
>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]
函数`flatten`可把一个矩阵转换为一个向量。
>flatten m;
[1,2,3,4,5,6]
检查维度 (shape)，行数(rows)和列数 (cols)：
>m=1..10$2:5;
>shape(m);
2 : 5
>rows(m)
2
>cols(m)
5
有两种获得某个单元格的值的方式：m.cell(row, col)和m[row, col]
>m=1..12$4:3;
>m;
#0 #1 #2
-- -- --
1 5 9
2 6 10
3 7 11
4 8 12
>m[1,2];
10
>m.cell(1,2);
10
有两种取得某一列的值的方式：一是m.col(index) ，其中index可以是数据对或者标量；二是m[index] 或 m[, index]，其中index可以是标量、数据对或者是向量。
>m=1..12$4:3;
>m;
#0 #1 #2
-- -- --
1 5 9
2 6 10
3 7 11
4 8 12
// 选出第1列生成一个向量。
>m[1];
[5,6,7,8]
// 选出第1列生成一个子矩阵。
>m[,1];
#0
--
5
6
7
8
// 选出第2列。
>m.col(2);
[9,10,11,12]
//若a<b，m[a:b]选出[a,b)范围内的列。
>m[1:3];
#0 #1
-- --
5 9
6 10
7 11
8 12
//若a>b，按倒序选出[b,a)范围内的列。
>m[2:0];
#0 #1
-- --
5 1
6 2
7 3
8 4
//单独选出某几列，此处为选出第0列和第2列。
>m[0 2];
#0 #1
-- --
1 9
2 10
3 11
4 12
有两种取得某一行的值的方式：一是m.row(index) ，其中index可以是数据对或者标量；二是m[index, ]，其中index可以是标量、数据对或者是向量。
>m=1..12$3:4;
>m;
#0 #1 #2 #3
-- -- -- --
1 4 7 10
2 5 8 11
3 6 9 12
// 返回第0行，类型为矩阵。
>m[0,];
#0 #1 #2 #3
-- -- -- --
1 4 7 10
// 用flatten函数将一个矩阵转化为向量。
>flatten(m[0,]);
[1,4,7,10]
// 选出第2行。
>m.row(2);
[3,6,9,12]
//若a<b，m[a:b,]选出[a,b)范围内的行。
>m[1:3, ];
#0 #1 #2 #3
-- -- -- --
2 5 8 11
3 6 9 12
//若a>b，按倒序选出[b,a)范围内的行。
>m[3:1, ];
#0 #1 #2 #3
-- -- -- --
3 6 9 12
2 5 8 11
///单独选出某几行，连续选出某些列。
>m[0 2,0..3];
#0 #1 #2 #3
-- -- -- --
1 4 7 10
3 6 9 12
有两种取得子矩阵的方式：第一种是m.slice(rowIndex,colIndex)或m[rowIndex,colIndex]，其中collndex和rowIndex是标量、向量或数据对。如果colIndex和rowIndex，那么上限边界值不包含在内，第二种与上述取行和取列方法类似。
>m=1..12$3:4;
>m;
#0 #1 #2 #3
-- -- -- --
1 4 7 10
2 5 8 11
3 6 9 12
// 选出第0,1行，第1,2列。
>m.slice(0:2,1:3);
#0 #1
-- --
4 7
5 8
// 选出[1,3)范围的行，[0,2)范围的列。
>m[1:3,0:2];
#0 #1
-- --
2 5
3 6
若追加一个向量到矩阵，该向量的长度必须是矩阵行数的倍数。
>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.
使用m[index]=X来修改某一列，X是一个标量或者向量。
>t1=1..50$10:5;
>t1;
#0 #1 #2 #3 #4
-- -- -- -- --
1 11 21 31 41
2 12 22 32 42
3 13 23 33 43
4 14 24 34 44
5 15 25 35 45
6 16 26 36 46
7 17 27 37 47
8 18 28 38 48
9 19 29 39 49
10 20 30 40 50
// 给第1列赋值200。
>t1[1]=200;
>t1;
#0 #1 #2 #3 #4
-- -- -- -- --
1 200 21 31 41
2 200 22 32 42
3 200 23 33 43
...
// 给第1列加200。
>t1[1]+=200;
>t1;
#0 #1 #2 #3 #4
-- -- -- -- --
1 400 21 31 41
2 400 22 32 42
3 400 23 33 43
...
// 将31..40的序列追加到第1列。
>t1[1]=31..40;
>t1;
#0 #1 #2 #3 #4
-- -- -- -- --
1 31 21 31 41
2 32 22 32 42
3 33 23 33 43
...
使用m[start:end] = X来修改多个列，X是一个标量或者向量。.
>t1=1..50$10:5;
>t1;
#0 #1 #2 #3 #4
-- -- -- -- --
1 11 21 31 41
2 12 22 32 42
3 13 23 33 43
4 14 24 34 44
5 15 25 35 45
6 16 26 36 46
7 17 27 37 47
8 18 28 38 48
9 19 29 39 49
10 20 30 40 50
// 追加向量101..130到第1,2,3列。
>t1[1:4]=101..130;
>t1;
#0 #1 #2 #3 #4
-- -- -- -- --
1 101 111 121 41
2 102 112 122 42
3 103 113 123 43
...
// 追加向量101..130到第3,2,1列。
>t1[4:1]=101..130;
>t1;
#0 #1 #2 #3 #4
-- -- -- -- --
1 121 111 101 41
2 122 112 102 42
3 123 113 103 43
...
// 第3,2,1列的每个元素加100。
>t1[4:1]+=100;
>t1;
#0 #1 #2 #3 #4
-- -- -- -- --
1 221 211 201 41
2 222 212 202 42
3 223 213 203 43
...
使用m[index,] = X修改某行，X是一个标量或者向量。
>t1=1..50$10:5;
// 赋值100到第1行。
>t1[1,]=100;
>t1;
#0 #1 #2 #3 #4
-- -- -- -- --
1 11 21 31 41
100 100 100 100 100
3 13 23 33 43
...
// 给第1行加100.
>t1[1,]+=100;
>t1;
#0 #1 #2 #3 #4
-- -- -- -- --
1 11 21 31 41
200 200 200 200 200
3 13 23 33 43
...
使用m[start:end,]=X修改多行，X是一个标量或者向量。
>t1=1..50$10:5;
// 将序列101..115赋值给第1,2,3行。
>t1[1:4,]=101..115;
>t1;
#0 #1 #2 #3 #4
--- --- --- --- ---
1 11 21 31 41
101 104 107 110 113
102 105 108 111 114
103 106 109 112 115
...
// 将序列101..115赋值给第3,2,1行。
>t1[4:1, ]=101..115;
>t1;
#0 #1 #2 #3 #4
--- --- --- --- ---
1 11 21 31 41
103 106 109 112 115
102 105 108 111 114
101 104 107 110 113
...
使用m[r1:r2, c1:c2]=X来修改矩阵窗口，X是一个标量或者向量。
>t1=1..50$5:10;
// 将序列101..110赋值给矩阵的第1,2行和第5~9列的窗口。
>t1[1:3,5:10]=101..110;
>t1;
#0 #1 #2 #3 #4 #5 #6 #7 #8 #9
-- -- -- -- -- --- --- --- --- ---
1 6 11 16 21 26 31 36 41 46
2 7 12 17 22 101 103 105 107 109
3 8 13 18 23 102 104 106 108 110
4 9 14 19 24 29 34 39 44 49
5 10 15 20 25 30 35 40 45 50
// 将序列101..110赋值给矩阵的第5~9行和第2~1列的窗口。
>t1=1..50$10:5;
>t1[5:10, 3:1]=101..110;
>t1;
#0 #1 #2 #3 #4
-- --- --- -- --
...
6 106 101 36 46
7 107 102 37 47
8 108 103 38 48
9 109 104 39 49
10 110 105 40 50
// 给矩阵的第9~5行，第1~2列加10。
>t1[10:5, 1:3]+=10;
>t1;
#0 #1 #2 #3 #4
-- --- --- -- --
1 11 21 31 41
2 12 22 32 42
3 13 23 33 43
4 14 24 34 44
5 15 25 35 45
6 116 111 36 46
7 117 112 37 47
8 118 113 38 48
9 119 114 39 49
10 120 115 40 50
可以使用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
DolphinDB中矩阵可被视为一个特殊的向量，所以大部分向量函数都适用于矩阵，比如求余弦值（cos）。当向量函数应用到矩阵时，结果仍然是一个矩阵。
>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中矩阵可以与其它矩阵或标量进行加减乘除运算。两个矩阵进行的运算为对应元素分别进行计算，所以两个矩阵的维度必须一致。
矩阵和标量进行运算:
>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的性能在一个数量级。
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
除了对矩阵进行所有元素进行计算外，DolphinDB支持对矩阵按行计算和按列计算。
支持按行计算的函数有:
- rowSum：计算矩阵每一行的和。
- rowSum2：计算矩阵每一行的平方和。
- rowVar：计算矩阵每一行的样本方差。
- rowMax：计算矩阵每一行的最大值。
- rowMin：计算矩阵每一行的最小值。
- rowProd：计算矩阵每一行元素的乘积。
- rowCount：计算矩阵每一行非NULL元素的个数。
- rowAvg：计算矩阵每一行的平均值。
- rowStd：计算矩阵每一行的样本标准差。
- rowXor：对矩阵每一行进行异或运算，若奇数列为true值，返回1，否则返回，计算中NULL值会被忽略。
- rowAnd：对矩阵每一行进行逻辑与的运算，计算中NULL值会被忽略。
- rowOr：对矩阵每一行进行或运算，若所有输入参数的所有列中至少有一列为true值，返回1，否则返回0，计算中NULL值会被忽略。
//计算矩阵每一行的和
>m=matrix([4.5 2.6 1.5, 1.5 4.8 5.9, 4.9 2.0 NULL])
rowSum(m);
[10.9,9.4,7.4]
采用高阶函数each对矩阵按列计算：
//计算矩阵每一列的和
>m=matrix([4.5 2.6 1.5, 1.5 4.8 5.9, 4.9 2.0 NULL])
>each(sum, m)
[8.6,12.2,6.9]
//两个矩阵按列进行矩阵相乘运算
>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
DolphinDB实现了以下矩阵分解函数：
lu(X,[permute=false])：对矩阵进行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
qr(X,[mode='full',pivoting=false])：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
函数svd(X, [fullMatrices=true,computeUV=true]计算矩阵 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]
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
schur(X, [sort]): 对矩阵进行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
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
函数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
pca(ds, [colNames], [k], [normalize], [maxIter]): 对数据源中指定列中的数据求主成分分析。
- ds是一个或多个数据源，通常由
sqlDS函数生成。 - colNames是字符串向量，表示数据源中的列名。默认值是数据源中所有列的列名。
- k是一个正整数，表示需要计算的主成分的个数。默认值是数据源中的列数。
- normalize是一个布尔值，表示是否除以标准差。默认值为false。
- maxIter是一个正整数，表示最大迭代次数。默认值为min(3*k,300)。
返回的结果是一个字典，包含以下键：
- 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]
从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函数中使用。