数据科学与机器学习（第 03 部分）：矩阵回归·进阶篇
（2/3）·从设计矩阵到 xTx 逆运算，解开多元动态回归卡了许久的编程死结
接上篇，我们继续深挖矩阵在回归里的实际用法。很多人卡在多元建模，是因为没把输入整理成设计矩阵，导致后续求逆和乘法全对不上号。先把第一列填 1 这件小事做对，后面几百个自变量才能动态塞进模型。
「矩阵积里的越界打印与极值钳制」
这段矩阵乘法内核里，作者用注释把两处调试打印藏了起来：一处盯 index 越界（index out），一处核对每一步累加结果。真在 MT5 里跑大规模矩阵时，这类 Print 一旦取消注释，日志会被刷爆，建议只在校验小规模样例时打开。
核心累加行 MultPl_Mat[index] += A[mat1_index] * B[mat2_index]; 之后紧跟 DBL_MAX_MIN(MultPl_Mat[index]);，作用是对单点结果做双精度上下界钳制，避免溢出变成 INF 或 denormal 拖慢后续求解。
收尾用 ArrayCopy(output_arr,MultPl_Mat); 把结果搬进输出数组，再 ArrayFree 释放临时矩阵——不 free 的话多次调用会吃满内存。下面贴出的 xTx 样例是 2×2 输出，对角 744.00000 与非对角 3257845.70000 可作为你本地验证矩阵转置乘自身的基准数值。
class=class="str">"cmt">//Print("index out ",index," index a ",mat1_index," index b ",mat2_index); MultPl_Mat[index] += A[mat1_index] * B[mat2_index]; DBL_MAX_MIN(MultPl_Mat[index]); class=class="str">"cmt">//Print(index," ",MultPl_Mat[index]); } ArrayCopy(output_arr,MultPl_Mat); ArrayFree(MultPl_Mat); } } k + (i*row2); j + (k*col2); Print("xTx"); MatrixPrint(xTx,tr_cols,tr_cols,class="num">5); xTx [ class="num">744.00000 class="num">3257845.70000 class="num">3257845.70000 class="num">14275586746.32998 ]
◍ 2x2 矩阵求逆的手动拆解与代码落地
做最小二乘线性回归时，xTx 这一步若卡在矩阵求逆上，手工算反而比调库更不容易出错。对 2x2 矩阵，逆的求法很直接：主对角线首尾互换，副对角线两个元素取负，再整体除以行列式 det = 主对角乘积 - 副对角乘积。 上面那段 MQL5 把这套动作写死了：只接 4 元素数组（即 2x2），超了就打印报错。实际跑出来，交换后的中间矩阵是 [14275586746, -3257846; -3257846, 744]，行列式算得 7477934261.0234375，最终逆矩阵打印为 [1.9090281, -0.0004357; -0.0004357, 0.0000001]。 外汇与贵金属行情噪声大，这类回归系数对样本区间极敏感，用逆矩阵算出的斜率只反映历史窗口内的线性倾向，实盘前务必在 MT5 用不同品种周期复算验证。 别把正态当圣经 代码里 DBL_MAX_MIN 那行是在给逆矩阵元素做溢出保护，复制去跑时若你自己的矩阵数量级差更大，det 接近 0 会直接把结果顶到极限值，这时候回归权重已经不可信。
class="type">void CSimpleMatLinearRegression::MatrixInverse(class="type">class="kw">double &Matrix[],class="type">class="kw">double &output_mat[]) { class=class="str">"cmt">// According to Matrix Rules the Inverse of a matrix can only be found when the class=class="str">"cmt">// Matrix is Identical Starting from a 2x2 matrix so this is our starting point class="type">int matrix_size = ArraySize(Matrix); if (matrix_size > class="num">4) Print("Matrix allowed using this method is a 2x2 matrix Only"); if (matrix_size==class="num">4) { MatrixtypeSquare(matrix_size); class=class="str">"cmt">//first step is we swap the first and the last value of the matrix class=class="str">"cmt">//so far we know that the last value is equal to arraysize minus one class="type">int last_mat = matrix_size-class="num">1; ArrayCopy(output_mat,Matrix); class=class="str">"cmt">// first diagonal output_mat[class="num">0] = Matrix[last_mat]; class=class="str">"cmt">//swap first array with last one output_mat[last_mat] = Matrix[class="num">0]; class=class="str">"cmt">//swap the last array with the first one class="type">class="kw">double first_diagonal = output_mat[class="num">0]*output_mat[last_mat]; class=class="str">"cmt">// second diagonal //adiing negative signs >>> output_mat[class="num">1] = - Matrix[class="num">1]; output_mat[class="num">2] = - Matrix[class="num">2]; class="type">class="kw">double second_diagonal = output_mat[class="num">1]*output_mat[class="num">2]; if (m_debug) { Print("Diagonal already Swapped Matrix"); MatrixPrint(output_mat,class="num">2,class="num">2); } class=class="str">"cmt">//formula for inverse is class="num">1/det(xTx) * (xtx)-class="num">1 class=class="str">"cmt">//determinant equals the product of the first diagonal minus the product of the second diagonal class="type">class="kw">double det = first_diagonal-second_diagonal; if (m_debug) Print("determinant =",det); for (class="type">int i=class="num">0; i<matrix_size; i++) { output_mat[i] = output_mat[i]*(class="num">1/det); DBL_MAX_MIN(output_mat[i]); } } } Diagonal already Swapped Matrix [ class="num">14275586746 -class="num">3257846 -class="num">3257846 class="num">744 ] determinant =class="num">7477934261.0234375 Print("inverse xtx"); MatrixPrint(inverse_xTx,class="num">2,class="num">2,_digits); class=class="str">"cmt">//inverse of simple lr will always be a 2x2 matrix
叉乘 xT 与 y 解出回归系数
把转置后的设计矩阵 xT 和因变量列 y 做矩阵乘法，得到 xTy 向量。实跑出来是 [10550016.7000000, 46241904488.2699585]，这两个数是后续求权重的中间量，不是价格，别拿去当信号。 用前面算好的 xTx 逆矩阵左乘 xTy，得到 Betas 系数向量，打印出来是 [-5524.40278, 4.49996]。第一位 -5524.40278 是常数项（Y 轴截距），第二位 4.49996 是斜率；这和第 01 部分用标量公式手算的结果一致，说明矩阵通路没接错。 常数项能落在 Betas[0]，是因为设计矩阵第一列全填了 1。这一列不是多余动作，它专门给截距留位，缺了这列逆矩阵维度都对不上。外汇与贵金属杠杆高、滑点乱，回归系数只是历史拟合，实盘照搬大概率失真。 下面这段代码就是上面两步的 MT5 落地，直接丢进带 MatrixMultiply / MatrixPrint 辅助函数的 EA 里能复现。
class="type">class="kw">double xTy[]; MatrixMultiply(xT,m_yvalues,xTy,tr_cols,tr_rows,tr_rows,class="num">1); class=class="str">"cmt">//class="num">1 at the end is because the y values matrix will always have one column which is it Print("xTy"); MatrixPrint(xTy,tr_rows,class="num">1,_digits); class=class="str">"cmt">//remember again??? how we find the output of our matrix row1 x column2 MatrixMultiply(inverse_xTx,xTy,Betas,class="num">2,class="num">2,class="num">2,class="num">1); class=class="str">"cmt">//inverse is a square 2x2 matrix while xty is a 2x1 Print("coefficients"); MatrixPrint(Betas,class="num">2,class="num">1,class="num">5); class=class="str">"cmt">// for simple lr our betas matrix will be a 2x1
「把多元回归塞进一个类里」
矩阵建模最实在的好处，是扩展维度时几乎不用动主逻辑，改几个参数就能从单变量切到多变量。难点一直卡在矩阵求逆，这块后面单说，眼下先把多元回归需要的骨架在独立函数库里搭出来。 我没把简单和多元写进同一个文件，而是拆开。只要你看懂了简单线性回归里那套计算，多元部分只是把 x 换成了设计矩阵，过程几乎一致。 Init() 里先抓 TestScript 选定的自变量列，塞进全局数组 m_XColsArray。用数组存列名，后面按正确顺序读 x 值会顺手很多。 数据集的每一行必须等长，只要有一行或一列错位，后面所有矩阵运算都会直接崩。随后把选中的 x 列拼成设计矩阵，因变量单独存进自己的矩阵。 设计矩阵第一列要填 1，这是截距项的位置。初始化时我把未转置的矩阵打印出来，方便你核对维度对不对。 下面这段是类声明的核心字段，跑通 multipleMapTregTestScript.mq5 后就能看到结构概貌：
class CMultipleMatLinearReg { class="kw">private: class="type">int m_handle; class="type">class="kw">string m_filename; class="type">class="kw">string DataColumnNames[]; class=class="str">"cmt">//store the column names from csv file class="type">int rows_total; class="type">int x_columns_chosen; class=class="str">"cmt">//Number of x columns chosen class="type">bool m_debug; class="type">class="kw">double m_yvalues[]; class=class="str">"cmt">//y values or dependent values matrix class="type">class="kw">double m_allxvalues[]; class=class="str">"cmt">//All x values design matrix class="type">class="kw">string m_XColsArray[]; class=class="str">"cmt">//store the x columns chosen on the Init class="type">class="kw">string m_delimiter; class="type">class="kw">double Betas[]; class=class="str">"cmt">//Array for storing the coefficients class="kw">protected: class="type">bool fileopen(); class="type">void GetAllDataToArray(class="type">class="kw">double& array[]); class="type">void GetColumnDatatoArray(class="type">int from_column_number, class="type">class="kw">double &toArr[]); class="kw">public: CMultipleMatLinearReg(class="type">void); ~CMultipleMatLinearReg(class="type">void);
◍ 多元线性回归的设计矩阵初始化细节
CMultipleMatLinearReg::Init 负责把外部 CSV 或缓存里的数据整理成可求逆的设计矩阵。它先把文件名、分隔符、调试开关存进成员变量，再用 StringSplit 按分隔符把 x_columns 字符串拆成列名数组，x_columns_chosen 直接等于该数组长度。 如果开了 debugmode，终端会打印「Init, number of X columns chosen =」并 ArrayPrint 出列名数组，方便你确认喂进去的因子数是否符合预期。 数据规整有个硬校验：用 rows_total 对 x_columns_chosen 取模，若余数不为 0 就 Alert 提示列长不一致，否则可能算出误导性系数。通过校验后，代码会造一个全 1 的 Temp_x 段拼到 m_allxvalues 头部，充当回归截距项。 MatrixUnTranspose 按 tr_cols = x_columns_chosen+1、tr_rows = single_rowsize 把扁平数组还原成矩阵，并把转置前副本存进 xT。外汇与贵金属数据常有缺口，跑这套前务必人工核对 rows_total 能被因子数整除，杠杆品种误算参数会引发实盘高风险。
class="type">void CMultipleMatLinearReg::Init(class="type">int y_column,class="type">class="kw">string x_columns="",class="type">class="kw">string filename=NULL,class="type">class="kw">string delimiter=",",class="type">bool debugmode=true) { class=class="str">"cmt">//--- pass some inputs to the global inputs since they are reusable m_filename = filename; m_debug = debugmode; m_delimiter = delimiter; class=class="str">"cmt">//--- class="type">class="kw">ushort separator = StringGetCharacter(m_delimiter,class="num">0); StringSplit(x_columns,separator,m_XColsArray); x_columns_chosen = ArraySize(m_XColsArray); ArrayResize(DataColumnNames,x_columns_chosen); class=class="str">"cmt">//--- if (m_debug) { Print("Init, number of X columns chosen =",x_columns_chosen); ArrayPrint(m_XColsArray); } class=class="str">"cmt">//--- GetAllDataToArray(m_allxvalues); GetColumnDatatoArray(y_column,m_yvalues); class=class="str">"cmt">// check for variance in the data set by dividing the rows total size by the number of x columns selected, there shouldn&class="macro">#x27;t be a reminder if (rows_total % x_columns_chosen != class="num">0) Alert("There are variance(s) in your dataset columns sizes, This may Lead to Incorrect calculations"); else { class=class="str">"cmt">//--- Refill the first row of a design matrix with the values of class="num">1 class="type">int single_rowsize = rows_total/x_columns_chosen; class="type">class="kw">double Temp_x[]; class=class="str">"cmt">//Temporary x array ArrayResize(Temp_x,single_rowsize); ArrayFill(Temp_x,class="num">0,single_rowsize,class="num">1); ArrayCopy(Temp_x,m_allxvalues,single_rowsize,class="num">0,WHOLE_ARRAY); class=class="str">"cmt">//after filling the values of one fill the remaining space with values of x class=class="str">"cmt">//Print("Temp x arr size =",ArraySize(Temp_x)); ArrayCopy(m_allxvalues,Temp_x); ArrayFree(Temp_x); class=class="str">"cmt">//we no longer need this array class="type">int tr_cols = x_columns_chosen+class="num">1, tr_rows = single_rowsize; ArrayCopy(xT,m_allxvalues); class=class="str">"cmt">//store the transposed values to their global array before we untranspose them MatrixUnTranspose(m_allxvalues,tr_cols,tr_rows); class=class="str">"cmt">//we add one to leave the space for the values of one if (m_debug) { Print("Design matrix"); MatrixPrint(m_allxvalues,tr_cols,tr_rows); } } }
设计矩阵回填与回归脚本落地
把 y 列抽进数组后，真正的坑在特征矩阵的构造：第一列必须全填 1，用来承载线性回归的截距项。代码里先用 single_rowsize = rows_total/x_columns_chosen 算每个特征列该摊多少行，再开临时数组 Temp_x 灌满 1，随后把原始 x 值拷贝进剩余位置，最后反置回 m_allxvalues 并释放临时内存。 MatrixUnTranspose(m_allxvalues,tr_cols,tr_rows) 这一步是核心，tr_cols 比选中的 x 列多 1，正是留给那列常数 1。若开了 m_debug，会直接 Print 出设计矩阵，肉眼核对首列是否全 1 比事后查系数更省时间。 实际跑一遍 OnStart，读入 NASDAQ_DATA.csv 并选 "1,3,4" 三列做特征，日志显示总数据量 2232 个、占 52 字节内存，设计矩阵首行是 [1, 4174, 13387, 35]、末段出现 [1, 4405, 14224, 56]。这种结构在 MT5 里验证一次，就能确认多维线性回归的输入没错位。 外汇与贵金属行情用同类矩阵建模时，样本外漂移倾向明显，高杠杆下误判概率会被放大，实盘前务必用历史分段回测。
GetColumnDatatoArray(y_column,m_yvalues); { class=class="str">"cmt">//--- Refill the first row of a design matrix with the values of class="num">1 class="type">int single_rowsize = rows_total/x_columns_chosen; class="type">class="kw">double Temp_x[]; class=class="str">"cmt">//Temporary x array ArrayResize(Temp_x,single_rowsize); ArrayFill(Temp_x,class="num">0,single_rowsize,class="num">1); ArrayCopy(Temp_x,m_allxvalues,single_rowsize,class="num">0,WHOLE_ARRAY); class=class="str">"cmt">//after filling the values of one fill the remaining space with values of x class=class="str">"cmt">//Print("Temp x arr size =",ArraySize(Temp_x)); ArrayCopy(m_allxvalues,Temp_x); ArrayFree(Temp_x); class=class="str">"cmt">//we no longer need this array class="type">int tr_cols = x_columns_chosen+class="num">1, tr_rows = single_rowsize; MatrixUnTranspose(m_allxvalues,tr_cols,tr_rows); class=class="str">"cmt">//we add one to leave the space for the values of one if (m_debug) { Print("Design matrix"); MatrixPrint(m_allxvalues,tr_cols,tr_rows); } } class="macro">#include "multipleMatLinearReg.mqh"; CMultipleMatLinearReg matreg; class=class="str">"cmt">//+------------------------------------------------------------------+ class=class="str">"cmt">//| Script program start function | class=class="str">"cmt">//+------------------------------------------------------------------+ class="type">void OnStart() { class=class="str">"cmt">//--- class="type">class="kw">string filename= "NASDAQ_DATA.csv"; matreg.Init(class="num">2,"class="num">1,class="num">3,class="num">4",filename); }