矩阵分解：更实用的建模

概述

大家好，欢迎阅读我的新文章，这篇文章再次以讲述为特色。

在上一篇文章"矩阵分解 ：基础知识"中 ，我谈到了读者们如何在一般计算中使用矩阵。不过，当时我想让大家明白计算是如何进行的，所以并没有太在意如何创建正确的矩阵模型。

您可能没有注意到，矩阵建模有点奇怪，因为只指定了列，而不是行和列。在阅读执行矩阵分解的代码时，这看起来非常奇怪。如果您希望看到列出的行和列，那么在尝试分解时可能会感到困惑。

此外，这种矩阵建模方法并不是最好的。这是因为当我们以这种方式对矩阵建模时，会遇到一些限制，迫使我们使用其他方法或函数，而如果以更合适的方式建模，这些方法或函数是不必要的。

虽然过程并不复杂，但要正确使用它，必须对其有充分的了解，因此我选择不在上一篇文章中详细介绍。在本文中，我们将更加冷静、不急不躁地考虑一切，以避免在矩阵建模方面产生错误的想法，从而确保矩阵的正确分解。

从本质上讲，矩阵是执行某些类型计算的更好方法，因为它们只需要我们做最少的额外执行工作。不过，出于同样的原因，在使用矩阵实现任何功能时都要非常小心。与传统建模不同的是，如果我们在矩阵建模中出现错误，就会得到相当奇怪的结果，顶多是在维护开发的代码时遇到很大困难。

矩阵建模时的注意事项

你可能会问，为什么我说上一篇文章中的建模很奇怪，因为这一切都说得通啊。也许，你理解了那篇文章中描述的一切。让我们看看上一篇文章中对矩阵进行建模的代码片段。也许这能让你更清楚地明白我要向你展示什么。下面是这段代码。

20. //+------------------------------------------------------------------+
21. void Arrow(const int x, const int y, const ushort angle, const uchar size = 100)
22. {
23.     double M_1[2][2] = {
24.                            cos(_ToRadians(angle)), sin(_ToRadians(angle)),
25.                           -sin(_ToRadians(angle)), cos(_ToRadians(angle))
26.                         },
27.            M_2[][2]  =  {
28.                            0.0,  0.0,
29.                            1.5, -.75,
30.                            1.0,  0.0,
31.                            1.5,  .75
32.                         },
33.            M_3[M_2.Size() / 2][2];

虽然它有效，但正如您从片段中看到的，这个问题非常令人困惑。这很难理解。在代码中编写多维数组时，我们使用以下配置：

Label[Dimension_01][Dimension_02]....[Dimension_N];

虽然这些代码在一般情况下看似简单，但当应用到矩阵分解领域时，这类符号就很难理解了。这是因为在处理矩阵时，这种在多维数组中的表现方式会造成混乱。你很快就会明白其中的原因。为简化起见，我们将把这一概念简化为二维系统，即只有行和列。不过，矩阵可以有任意维数。如果简化为二维矩阵，其符号将如下所示：

这是何意？字母 "M" 代表矩阵名称，"l" 代表行数，"c" 代表列数。一般来说，在提及矩阵元素时，我们会使用相同的符号，如下所示：

请注意，矩阵 M 现在由六个元素组成，可以使用行和列系统进行访问。但请仔细观察上一张图片：它看起来像什么？根据您的经验，它可能看起来有所不同。然而，用编程术语来说，这是什么意思呢？所以说，这是一个二维数组。但采用这种方法后，矩阵看起来是静态的。而矩阵是一个动态实体，可以根据我们计算的内容以不同的方式读取。在某些情况下，我们需要对角阅读，在其他情况下，我们需要逐列阅读，也需要逐行阅读。因此，将矩阵视为多维数组会将动态的东西变成静态的东西。这使得某些矩阵计算的实现变得复杂。

不过，我们不是来学习编程基础知识的。亲爱的读者们，我想告诉你们，如何以最佳方式将任何矩阵可视化，同时让它保持其应有的动态实体。

因此，将矩阵视为数组的想法并没有错。当把它看作一个多维数组时，问题就出现了。理想情况下，应将矩阵视为一维数组，但没有固定的宽度。这一点听起来可能很奇怪，但问题在于语言。目前还没有合适的术语或表达方式来描述某些概念。

还记得艾萨克-牛顿第一次解释他的计算时所面临的困难吗？由于某些概念缺乏适当的术语，这肯定相当困难。不过，还是回到我们的问题上来吧。

我们已经有了将矩阵表示为数组的方法。但现在出现了第一个问题：如何对数组中的元素排序？如果你不明白这为什么是个问题，那么你就不明白有关矩阵分解的其他相关问题。元素的排序方式在很大程度上决定了分解代码的实现方式。虽然两种情况下的结果是一样的，但我们可能需要在数组内进行更多或更少的转换才能得到正确的结果。这是因为计算机内存是线性的。如果我们将元素排序为 t11、t12、t21、t22、t31 和 t32（这在很多情况下都没问题），那么我们就需要强制跳转，以便以不同的顺序访问元素来执行特定计算。无论我们做什么，这种情况肯定会发生。

通常，当我们对数组元素进行排序时，是通过逐行读取来实现的。但是，您可以以不同的方式对它们进行排序：例如，逐列阅读或以任何其他方式自行决定。你可以随心所欲地去做。

一旦这样做了，就会出现另一个小问题，这个问题与编程无关，而是与如何进行矩阵分解有关。在这种情况下，我就不细说了。原因很简单：如何处理这个问题，取决于我们如何实现矩阵或使用哪种矩阵。我建议学习相关书籍或数学文献，以了解如何解决可能出现的问题，并正确地进行分解。然而，实现代码从来都不是一件难事，因为代码是开发过程中最简单的部分。最难的部分是了解如何进行分解。上图中矩阵行列式的计算就是一个例子。请注意，这不是一个方阵，在这种情况下构造行列式与我们在方阵的情况下所做的略有不同。如前所述，有必要研究数学部分，因为每个具体情况都可能需要特殊的方法。

混乱

好了，现在一切都解释清楚了，我们可以进入编码阶段了。让我们继续上一篇文章中讨论的例子，因为这个例子非常实用，也很容易理解。此外，这个例子还很容易说明在某些情况下使用矩阵计算的另一个优势。让我们回忆一下代码的样子。下面是这段代码：

01. //+------------------------------------------------------------------+
02. #property copyright "Daniel Jose"
03. #property indicator_chart_window
04. #property indicator_plots 0
05. //+------------------------------------------------------------------+
06. #include <Canvas\Canvas.mqh>
07. //+------------------------------------------------------------------+
08. CCanvas canvas;
09. //+------------------------------------------------------------------+
10. #define _ToRadians(A) (A * (M_PI / 180.0))
11. //+------------------------------------------------------------------+
12. void MatrixA_x_MatrixB(const double &A[][], const double &B[][], double &R[][], const int nDim)
13. {
14.     for (int c = 0, size = (int)(B.Size() / nDim); c < size; c++)
15.     {
16.         R[c][0] = (A[0][0] * B[c][0]) + (A[0][1] * B[c][1]);
17.         R[c][1] = (A[1][0] * B[c][0]) + (A[1][1] * B[c][1]);
18.     }
19. }
20. //+------------------------------------------------------------------+
21. void Arrow(const int x, const int y, const ushort angle, const uchar size = 100)
22. {
23.     double M_1[2][2]{
24.                         cos(_ToRadians(angle)), sin(_ToRadians(angle)),
25.                         -sin(_ToRadians(angle)), cos(_ToRadians(angle))
26.                     },
27.            M_2[][2] {
28.                         0.0,  0.0,
29.                         1.5, -.75,
30.                         1.0,  0.0,
31.                         1.5,  .75
32.                     },
33.            M_3[M_2.Size() / 2][2];
34. 
35.     int dx[M_2.Size() / 2], dy[M_2.Size() / 2];
36.     
37.     MatrixA_x_MatrixB(M_1, M_2, M_3, 2);
38.     ZeroMemory(M_1);
39.     M_1[0][0] = M_1[1][1] = size;
40.     MatrixA_x_MatrixB(M_1, M_3, M_2, 2);
41. 
42.     for (int c = 0; c < (int)M_2.Size() / 2; c++)
43.     {
44.         dx[c] = x + (int) M_2[c][0];
45.         dy[c] = y + (int) M_2[c][1];
46.     }
47. 
48.     canvas.FillPolygon(dx, dy, ColorToARGB(clrPurple, 255));
49.     canvas.FillCircle(x, y, 5, ColorToARGB(clrRed, 255));
50. }
51. //+------------------------------------------------------------------+
52. int OnInit()
53. {    
54.     int px, py;
55.     
56.     px = (int)ChartGetInteger(0, CHART_WIDTH_IN_PIXELS, 0);
57.     py = (int)ChartGetInteger(0, CHART_HEIGHT_IN_PIXELS, 0);
58. 
59.     canvas.CreateBitmapLabel("BL", 0, 0, px, py, COLOR_FORMAT_ARGB_NORMALIZE);
60.     canvas.Erase(ColorToARGB(clrWhite, 255));
61.         
62.     Arrow(px / 2, py / 2, 160);
63. 
64.     canvas.Update(true);
65.     
66.     return INIT_SUCCEEDED;
67. }
68. //+------------------------------------------------------------------+
69. int OnCalculate(const int rates_total, const int prev_calculated, const int begin, const double &price[])
70. {
71.     return rates_total;
72. }
73. //+------------------------------------------------------------------+
74. void OnDeinit(const int reason)
75. {
76.     canvas.Destroy();
77. }
78. //+------------------------------------------------------------------+

由于此代码已在上一篇文章中描述，并且也在其附录中提供，因此我们在此不再重复。让我们只关注 Arrow 过程和 MatrixA_x_MatrixB 过程，它们实际上在我们的演示代码中做了一些事情。

首先，让我们改变创建数组的方式。这些变化实际上比看起来更简单、更微妙，但它们仍然会改变我们开始的方式。

让我们从上面代码中显示的 M_1 数组开始。研究它，了解它的内容。

请注意，这里有一些有趣的东西，这与上一节讨论的主题完全相关。第 23 行声明了一个两行两列的数组或矩阵 M_1。在这种情况下，这很简单。这同样适用于数组或矩阵 M_2，它没有行数定义，但有列数定义，列数等于 2。到目前为止，没有任何疑问或问题；很容易理解我们在做什么，因为代码中的每一行新行都对应于矩阵中的新行，而该行中的每一列都将是矩阵中的新列。这一切都很简单。然而，问题出在 MatrixA_x_MatrixB 过程的代码中。你能找出问题所在吗？不能？还是可以？也许，你可以看到它，但不明白。

事实上，问题出在负责矩阵乘法的 MatrixA_x_MatrixB 过程上。由于我们在这里做的事情非常简单，因此不会影响代码的实现。但是，我们仍然需要解决这个问题，因为它确实存在，而且阻碍了我们的进步。一旦我们需要执行更复杂的乘法运算，这就会成为一个真正的障碍。

如果你没有注意到这个问题，请注意变量 B。它是一个将另一个矩阵相乘的矩阵。为了更好地可视化问题，让我们将矩阵 B 缩减为一列但有两行，因为矩阵 A 有两列。要进行乘法运算，您需要做下图所示的操作。

这样就能得出正确的结果。这正是对矩阵进行乘法或分解时发生的情况。如需了解更详细的信息，我建议您研究一下这个问题。在这两种情况下，乘法运算都是在列中逐行进行的，结果是在列中保留数值。出现这种情况是因为第二个因子是一个列。如果乘法以列的形式进行，则所得值必须以行的形式放置。

虽然说了这么多似乎只会让你感到困惑，但事实上，这正是矩阵乘法的规则。当其他程序员分析代码时，他们不会发现所有操作都是正确的。这是第 16 行和第 17 行中可以看到的问题，其中必须执行上图所示的计算。但在矩阵 A 中，我们使用了逐行计算，而在矩阵 B 中，我们也使用了逐行计算。这让我们感到困惑，尽管在这个简单的例子中，它提供了正确的数据。

你可能会想：那又怎样？如果一切正常，就保持原样。但这不对。如果在某个时候你需要改进代码来处理更大的矩阵或覆盖稍微不同的情况，那么我们最终会试图让一些相对简单的东西工作，从而使问题变得更加复杂。因此，这些事情应该做得很好，以避免混淆。为了解决这个问题，我们需要用一个稍有不同的模型来取代这个模型。虽然一开始可能看起来很混乱，但正如上图所示，这是正确的。

为了避免混淆，让我们看看如何在新主题中做到这一点。

让事情更简单

虽然我们在这里要做的事情可能看起来很复杂，但如果做得正确和谨慎，它真的不是。这可能更乏味，但我们将以一种相当简单的方式完成。请看下图。

这种表示法表明这两个矩阵是等价的。这是因为左侧有元素 A 的矩阵是行矩阵。当将其转换为右侧有 B 个元素的矩阵时，在这种情况下是列矩阵，我们不必做太多更改，在某些情况下可能只是旋转或更改索引。然而，我们可能需要多次这样做：将行矩阵转换为列矩阵，反之亦然，对矩阵进行分解。但要使这个例子正确，a11 必须等于 b11，a21 必须等于 b12。

在某些情况下，您可能会注意到，无需进行任何重大更改，只需更改索引即可将行结构转换为列之一。但如何正确地做到这一点，要视具体情况而定。因此，从相关数学书籍中了解矩阵分解的方法非常重要。

现在让我们回到我们在上一个主题中看到的代码。首先，让我们纠正一下表述。让我们更改代码，如下所示。

20. //+------------------------------------------------------------------+
21. void Arrow(const int x, const int y, const ushort angle, const uchar size = 100)
22. {
23.     double  M_1[]= {
24.                        cos(_ToRadians(angle)), sin(_ToRadians(angle)),
25.                        -sin(_ToRadians(angle)), cos(_ToRadians(angle))
26.                    },
27.             M_2[]= {
28.                        0.0,  0.0,
29.                        1.5, -.75,
30.                        1.0,  0.0,
31.                        1.5,  .75
32.                    },
33.            M_3[M_2.Size()];

注意，数组中不再有索引。矩阵结构如代码所示创建，即行是行，列是列。但要注意 M_3 等动态数组。我们必须正确定义分配的长度，以避免在访问这些矩阵的元素时出现运行时错误。

这是第一部分。在此之后，过程 MatrixA_x_MatrixB 将不再能够理解如何进行计算。此外，第 42 行的循环（见上一主题中的代码）将无法正确解码数据。

首先，我们将修正在 MatrixA_x_MatrixB 过程中进行的计算，然后修正绘图。不过，在我们继续之前，请注意矩阵 M_1 几乎没有变化，因为它是对称的正方形。这是我们可以在行和列中读取它的情况之一。

由于大部分代码都会发生变化，让我们来看看这段新代码。我认为这将更容易解释这一切将如何运作。

09. //+------------------------------------------------------------------+
10. #define _ToRadians(A) (A * (M_PI / 180.0))
11. //+------------------------------------------------------------------+
12. void MatrixA_x_MatrixB(const double &A[], const ushort Rows, const double &B[], double &R[])
13. {
14.     uint Lines = (uint)(A.Size() / Rows);
15. 
16.     for (uint cbl = 0; cbl < B.Size(); cbl += Rows)
17.         for (uint car = 0; car < Rows; car++)
18.         {
19.             R[car + cbl] = 0;
20.             for (uint cal = 0; cal < Lines; cal++)
21.                 R[car + cbl] += (A[(cal * Rows) + car] * B[cal + cbl]);
22.         }
23. }
24. //+------------------------------------------------------------------+
25. void Plot(const int x, const int y, const double &A[])
26. {
27.     int dx[], dy[];
28. 
29.     for (uint c0 = 0, c1 = 0; c1 < A.Size(); c0++)
30.     {
31.         ArrayResize(dx, c0 + 1, A.Size());
32.         ArrayResize(dy, c0 + 1, A.Size());
33.         dx[c0] = x + (int)(A[c1++]);
34.         dy[c0] = y + (int)(A[c1++]);
35.     }
36. 
37.     canvas.FillPolygon(dx, dy, ColorToARGB(clrPurple, 255));
38.     canvas.FillCircle(x, y, 5, ColorToARGB(clrRed, 255));
39. 
40.     ArrayFree(dx);
41.     ArrayFree(dy);
42. }
43. //+------------------------------------------------------------------+
44. void Arrow(const int x, const int y, const ushort angle, const uchar size = 100)
45. {
46.     double  M_1[]={
47.                      cos(_ToRadians(angle)), sin(_ToRadians(angle)),
48.                      -sin(_ToRadians(angle)), cos(_ToRadians(angle))
49.                   },
50.             M_2[]={
51.                      0.0,  0.0,
52.                      1.5, -.75,
53.                      1.0,  0.0,
54.                      1.5,  .75
55.                   },
56.            M_3[M_2.Size()];
57.     
58.     MatrixA_x_MatrixB(M_1, 2, M_2, M_3);
59.     ZeroMemory(M_1);
60.     M_1[0] = M_1[3] = size;
61.     MatrixA_x_MatrixB(M_1, 2, M_3, M_2);
62.     Plot(x, y, M_2);
63. }
64. //+------------------------------------------------------------------+

现在一切似乎都变得更加复杂了。我们真的需要这种混乱吗？冷静下来，我亲爱的读者。该代码既不复杂也不混乱。它其实非常简单明了。此外，它还非常灵活。我会解释我为什么这么想。

在这个片段中，两个矩阵之间的乘法是正确完成的。在这种情况下，我们实现了一个行中列模型，从而得到一个新行作为乘法的乘积。乘法正好发生在第 21 行。我知道第 12 行执行的过程可能看起来很复杂，但事实并非如此。所有这些变量的存在都是为了确保正确访问数组中的特定索引。

由于没有测试来确保一个矩阵中的行数等于另一个矩阵的列数，因此在使用此过程时必须小心，以避免运行时错误。

现在，必须按列组织的矩阵 A 乘以必须按行组织的矩阵 B，并将结果放入按行排列的矩阵 R 中，我们可以看到调用是如何进行的。这发生在第 58 和 61 行。首先来看第 58 行。第一个要传递的矩阵是列矩阵 M_1。作为第二个矩阵，我们传递矩阵 M_2，它是动态的，按行排列。如果我们在第 58 行改变这个顺序，计算值就会出错。我知道，2 x 3 与 3 x 2 不同，这看起来很奇怪，但说到矩阵，因数的顺序确实会影响结果。但请注意，这与之前的做法类似。当然，调用的第二个参数指定了第一个矩阵中的列数。由于我们在这里的目的是说明，我将不检查矩阵 B（在本例中为 M_2）与矩阵 A（M_1）相乘是否可能。我假设你不会改变结构，也许只是增加矩阵 M_2 的行数，这不会影响 MatrixA_x_MatrixB 中的分解。

好，第 59 行和第 60 行已经在前一篇文章中解释过了，第 61 行发生的事情等同于第 58 行发生的事情，因此我们可以认为解释已经给出。不过，在第 62 行中，我们有一个新的调用，它调用了第 25 行中的过程。第 25 行可用来直接在屏幕上绘制其他对象或图像。我们需要做的就是给这个过程传递一个行类型矩阵，其中每一行代表要绘制的图像的一个顶点。矩阵中应包含的点应按以下顺序排列：X 在第一列，Y 在第二列，因为我们是在二维屏幕上绘图。因此，即使我们在三维中绘制东西，我们也需要在矩阵中进行一些计算，以便将 Z 维度中的顶点投影到 XY 平面上。我们不必担心分离问题或矩阵的大小。Plot 方法就可以做到这一点。请注意，在第 27 行中，我们声明了两个动态数组。在第 29 行，我们创建了一个有趣的 for 循环，它不同于许多人习惯的循环编写方式。在这个循环中，我们有两个变量：一个由循环控制，另一个在代码内部。变量 c0 由一个循环管理，该循环会尝试正确实例化 dx 和 dy 的索引，并在必要时分配更多内存。这是在第 31 和 32 行完成的。虽然这看起来效率不高，但看看 MQL5 文档，你就会发现情况并非如此。

另一方面，变量 c1 是在代码中管理的。为了理解这一点，让我们看看第 33 行和第 34 行中的 c1 发生了什么。每迭代一次就增加一。不过，决定 for 循环何时终止的是 c1。查看第 29 行的循环终止条件。

这里我们不检查数组或矩阵 A 的结构是否正确。如果一行中没有两列，那么在for循环的某个时刻，将生成运行时错误，表示有人试图访问超出范围的位置。因此，在创建矩阵 A 时应小心谨慎。

最后，我们使用第 37 行调用 CCanvas 类中的一个方法。这样做是为了在图表上绘制矩阵形状。第 38 行显示绘图起始点的位置。如果一切顺利，我们将看到下图。

当然，在第 40 行和第 41 行中，我们将已分配的内存返回给操作系统，从而释放了已分配的资源。

最后的探讨

在过去的两篇文章中，我试图介绍一些很多人都觉得困难的内容：矩阵分解。我知道，这里介绍的材料并没有涵盖在计算中使用矩阵的所有方面和好处。然而，我想强调的是，了解如何开发计算矩阵的方法非常重要。三维程序，包括游戏和矢量图像编辑器，都使用这种分解。即使是看似简单的程序（如光栅图像程序）也会使用矩阵计算来优化系统性能。

在这两篇文章中，我试图用最简单的方式来说明这一点，因为我注意到，许多编程新手根本就忽视了学习某些知识的必要性。他们经常使用库或其他花哨的软件，却不知道背后发生了什么。

在这里，我尽量简化流程。我只关注一个操作：两个矩阵相乘。为了展示它有多简单，我没有展示如何使用 GPU 来完成 MatrixA_x_MatrixB 过程中的工作，而这正是许多情况下会发生的事情。在这里，我们只使用 CPU。但一般来说，矩阵分解通常不在 CPU 上进行，而是在 GPU 上进行，使用 OpenCL 进行正确编程。

我认为，在某种程度上，我已经证明，我们可以做得比某些人认为的要多得多。研究文章中介绍的材料，因为很快我将讨论一个矩阵分解是基础的话题。

我想补充的最后一个细节是，两个矩阵相乘的目的是解决一个特定的问题。这与书本或学校传统上教授的方法相去甚远。一般来说，不要指望使用这种因子分解模型，因为它不起作用。它只适用于解决所提出的问题，即旋转物体。