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