LISP - 矩阵乘法

Lisp 项目：实现矩阵乘法

矩阵乘法是计算机图形学、科学计算和机器学习等领域中的一个基本操作。本章将指导你如何在 Common Lisp 中实现它，作为一个实用的学习练习。

免责声明：对于生产级的科学计算，强烈建议使用专门优化的库，如 Magicl，它可以通过 Quicklisp 安装。我们这里的实现仅用于教育目的。

数据表示：使用二维数组

尽管你可以使用列表的列表来表示矩阵，但在 Common Lisp 中，更自然、性能更好的选择是使用二维数组 (2D array)。数组提供了快速、常数时间的元素访问，这对于数值计算至关重要。

;; 创建一个 2x3 矩阵（2 行，3 列），用零初始化
(make-array '(2 3) :initial-element 0)

;; 访问第 1 行、第 2 列的元素
(let ((matrix (make-array '(2 3) :initial-contents '((1 2 3) (4 5 6)))))
  (aref matrix 1 2)) ; 返回 6

矩阵乘法的逻辑

要将矩阵 A（维度 m x n）乘以矩阵 B（维度 n x p）：

• 兼容性：矩阵 A 的列数 (n) 必须等于矩阵 B 的行数 (n)。
• 结果维度：结果矩阵 C 的维度将是 m x p。
• 计算：每个元素 C(i, j) 是矩阵 A 的第 i 行与矩阵 B 的第 j 列的点积。

实现

我们将编写一个 multiply-matrices 函数，它接受两个二维数组作为参数，并返回它们相乘后的新二维数组。我们将使用 loop 宏来实现一个清晰高效的乘法。

matrix-ops.lisp

(defun multiply-matrices (matrix-a matrix-b)
  "Multiplies matrix-a by matrix-b and returns the resulting matrix."
  (let* ((dims-a (array-dimensions matrix-a))
         (dims-b (array-dimensions matrix-b))
         (rows-a (first dims-a))
         (cols-a (second dims-a))
         (rows-b (first dims-b))
         (cols-b (second dims-b)))

    ;; Check for dimension compatibility
    (unless (= cols-a rows-b)
      (error "Matrix dimensions are incompatible for multiplication."))

    ;; Create the result matrix, initialized to 0.0 for floating point precision
    (let ((result (make-array (list rows-a cols-b) :initial-element 0.0)))
      (loop for i from 0 below rows-a
            do (loop for j from 0 below cols-b
                     do (loop for k from 0 below cols-a
                              do (incf (aref result i j)
                                     (* (aref matrix-a i k)
                                        (aref matrix-b k j))))))
      result)))

;; --- 示例用法 ---
(let ((matrix-a (make-array '(2 3) :initial-contents '((1 2 3)
                                                      (4 5 6))))
      (matrix-b (make-array '(3 2) :initial-contents '((7 8)
                                                      (9 10)
                                                      (11 12)))))
  (let ((result (multiply-matrices matrix-a matrix-b)))
    (format t "Matrix A:~%~a~%" matrix-a)
    (format t "Matrix B:~%~a~%" matrix-b)
    (format t "Result C:~%~a~%" result)))

输出

Matrix A:
#2A((1 2 3) (4 5 6))
Matrix B:
#2A((7 8) (9 10) (11 12))
Result C:
#2A((58.0 64.0) (139.0 154.0))

测试我们的函数

好的代码应该是可测试的。一个简单的测试方法是使用断言 (assertion)。例如，任何矩阵乘以单位矩阵都应该得到原始矩阵。我们可以使用 assert 来验证这一点。

(defun test-matrix-multiplication ()
  "A simple test case for our function."
  (let* ((matrix (make-array '(2 3) :initial-contents '((1 2 3) (4 5 6))))
         (identity-matrix (make-array '(3 3) :initial-contents '((1 0 0)
                                                                  (0 1 0)
                                                                  (0 0 1))))
         (result (multiply-matrices matrix identity-matrix)))
    ;; equalp 按元素比较数组。
    (assert (equalp result matrix))
    (format t "Test passed!~%")))

;; 运行测试
(test-matrix-multiplication)