From 9d7a7fee9ee095168ea2435d3e211a7b1f71f832 Mon Sep 17 00:00:00 2001 From: Brian Eberman Date: Fri, 29 Jul 2022 14:52:13 -0400 Subject: [PATCH 01/12] extended dot product to support dot product with a vector with itself. This was blocked by the original design --- lisp/data-vector.lisp | 12 ++++++++++++ lisp/sequence/vector.lisp | 1 + test/data-vector.lisp | 11 +++++++++++ 3 files changed, 24 insertions(+) diff --git a/lisp/data-vector.lisp b/lisp/data-vector.lisp index de80077..22395bb 100644 --- a/lisp/data-vector.lisp +++ b/lisp/data-vector.lisp @@ -434,3 +434,15 @@ applying the function to each element of the vectors." ((vector1 row-vector) (vector2 column-vector) &optional scalar) "Return the dot product of vector1 and vector2." (inner-product-vector (contents vector1) (contents vector2) scalar)) + +(defmethod product :before + ((vector1 data-vector) (vector2 data-vector) &optional scalar) + "Verify that the dimensions are equal." + (declare (ignore scalar)) + (unless (= (vector-length vector1) (vector-length vector2)) + (error "VECTOR1 and VECTOR2 are not of equal length."))) + +(defmethod product + ((vector1 data-vector) (vector2 data-vector) &optional scalar) + "Return the dot product of vector1 and vector2." + (inner-product-vector (contents vector1) (contents vector2) scalar)) diff --git a/lisp/sequence/vector.lisp b/lisp/sequence/vector.lisp index dace026..75619de 100644 --- a/lisp/sequence/vector.lisp +++ b/lisp/sequence/vector.lisp @@ -110,3 +110,4 @@ (inner-product-vector vector1 vector2 scalar) (error "VECTOR1(~D) and VECTOR2(~D) are not of equal length." (length vector1) (length vector2)))) + diff --git a/test/data-vector.lisp b/test/data-vector.lisp index e0b3589..9da481f 100644 --- a/test/data-vector.lisp +++ b/test/data-vector.lisp @@ -1222,3 +1222,14 @@ (linear-algebra:product (linear-algebra:column-vector 1 2 3) (linear-algebra:row-vector 1 2 3)))) + +;; test the inner product +(define-test inner-product-vector + (:tag :data-vector :product) + ;; Real vectors + (assert-rational-equal + 14 + (linear-algebra-kernel:inner-product-vector + (linear-algebra:row-vector 1 2 3) + (linear-algebra:row-vector 1 2 3))) +) From fe9c4d0509af514ceaa7cbfcf9203460023adbba Mon Sep 17 00:00:00 2001 From: Brian Eberman Date: Thu, 4 Aug 2022 11:49:46 -0400 Subject: [PATCH 02/12] added equality testing --- lisp/dense-matrix.lisp | 5 +++++ lisp/interface/fundamental-ops.lisp | 4 ++++ lisp/linear-algebra.lisp | 1 + 3 files changed, 10 insertions(+) diff --git a/lisp/dense-matrix.lisp b/lisp/dense-matrix.lisp index f06b28b..4d3d1c8 100644 --- a/lisp/dense-matrix.lisp +++ b/lisp/dense-matrix.lisp @@ -34,6 +34,7 @@ (:documentation "Dense matrix object.")) + (defmethod initialize-instance :after ((self dense-matrix) &rest initargs &key dimensions element-type initial-element initial-contents) @@ -110,6 +111,10 @@ "Return true if object is a dense matrix." (typep object 'dense-matrix)) +(defmethod mat-equal + ((matrix1 dense-matrix) (matrix2 dense-matrix)) + (equalp (contents matrix1) (contents matrix2))) + (defmethod matrix-in-bounds-p ((matrix dense-matrix) (row integer) (column integer)) "Return true if row and column do not exceed the dimensions of matrix." diff --git a/lisp/interface/fundamental-ops.lisp b/lisp/interface/fundamental-ops.lisp index 3d142d1..44aab03 100644 --- a/lisp/interface/fundamental-ops.lisp +++ b/lisp/interface/fundamental-ops.lisp @@ -92,3 +92,7 @@ (defgeneric ninvert (matrix) (:documentation "Return the invert of the matrix with in-place decomposition.")) + +(defgeneric mat-equal (matrix1 matrix2) + (:documentation + "Return if the two arrays elments are the same")) diff --git a/lisp/linear-algebra.lisp b/lisp/linear-algebra.lisp index 1586a88..4062d58 100644 --- a/lisp/linear-algebra.lisp +++ b/lisp/linear-algebra.lisp @@ -65,6 +65,7 @@ :matrix-column-dimension :mref :copy-matrix + :mat-equal :submatrix :replace-matrix :matrix-validated-range) From 40b3ba1e97373912de4331499bad116b6e21c598 Mon Sep 17 00:00:00 2001 From: Brian Eberman Date: Sun, 28 Aug 2022 10:15:48 -0400 Subject: [PATCH 03/12] adding functions for equality and pointwise operations --- lisp/data-vector.lisp | 58 ++++++++++++++++++++++++++ lisp/interface/fundamental-ops.lisp | 9 ++-- lisp/interface/matrix.lisp | 4 ++ lisp/interface/vector.lisp | 16 +++++++ lisp/kernel/binary-operations.lisp | 24 ++++++++++- lisp/kernel/linear-algebra-kernel.lisp | 2 + lisp/linear-algebra.lisp | 5 +++ 7 files changed, 114 insertions(+), 4 deletions(-) diff --git a/lisp/data-vector.lisp b/lisp/data-vector.lisp index 22395bb..97245c8 100644 --- a/lisp/data-vector.lisp +++ b/lisp/data-vector.lisp @@ -200,6 +200,9 @@ applying the function to each element of the vectors." function vectors)) +(defmethod reduce-vector ((function function) (vector data-vector) initial-value) + (reduce function (contents vector) :initial-value initial-value)) + ;;; Data vector transformations (defmethod apply-rotation :before @@ -446,3 +449,58 @@ applying the function to each element of the vectors." ((vector1 data-vector) (vector2 data-vector) &optional scalar) "Return the dot product of vector1 and vector2." (inner-product-vector (contents vector1) (contents vector2) scalar)) + +(defmethod vec-equal :before + ((vector1 data-vector) (vector2 data-vector)) + "Verify that the dimensions are equal." + (declare (ignore scalar)) + (unless (= (vector-length vector1) (vector-length vector2)) + (error "VECTOR1 and VECTOR2 are not of equal length."))) + +(defmethod vec-equal + ((vector1 data-vector) (vector2 data-vector)) + (equalp (contents vector1) (contents vector2))) + + +(defmethod elem-divide :before + ((vector1 data-vector) (vector2 data-vector)) + "Verify that the dimensions are equal." + (declare (ignore scalar)) + (unless (= (vector-length vector1) (vector-length vector2)) + (error "VECTOR1 and VECTOR2 are not of equal length."))) + +(defmethod elem-divide + ((vector1 data-vector) (vector2 data-vector)) + (make-instance + (common-class-of vector1 vector2) + :contents + (element-divide-vector + (contents vector1) (contents vector2)) + )) + +(defmethod elem-multiply :before + ((vector1 data-vector) (vector2 data-vector)) + "Verify that the dimensions are equal." + (declare (ignore scalar)) + (unless (= (vector-length vector1) (vector-length vector2)) + (error "VECTOR1 and VECTOR2 are not of equal length."))) + +(defmethod elem-multiply + ((vector1 data-vector) (vector2 data-vector)) + (make-instance + (common-class-of vector1 vector2) + :contents + (element-multiply-vector + (contents vector1) (contents vector2)) + )) + +(defmethod distance :before + ((vector1 data-vector) (vector2 data-vector) &optional (measure 2)) + "Verify that the dimensions are equal." + (declare (ignore scalar)) + (unless (= (vector-length vector1) (vector-length vector2)) + (error "VECTOR1 and VECTOR2 are not of equal length."))) + +(defmethod distance + ((vector1 data-vector) (vector2 data-vector) &optional (measure 2)) + (norm (subtract vector1 vector2) measure)) diff --git a/lisp/interface/fundamental-ops.lisp b/lisp/interface/fundamental-ops.lisp index 44aab03..f438f71 100644 --- a/lisp/interface/fundamental-ops.lisp +++ b/lisp/interface/fundamental-ops.lisp @@ -32,6 +32,10 @@ (:documentation "Return the norm according to measure.")) +(defgeneric distance (vector1 vector2 &optional measure) + (:documentation + "Return the distance between vector1 and vector2 using the measure")) + (defgeneric transpose (vector-or-matrix) (:documentation "Transpose the vector or matrix.")) @@ -93,6 +97,5 @@ (:documentation "Return the invert of the matrix with in-place decomposition.")) -(defgeneric mat-equal (matrix1 matrix2) - (:documentation - "Return if the two arrays elments are the same")) + + diff --git a/lisp/interface/matrix.lisp b/lisp/interface/matrix.lisp index a8f84fa..f01e9e7 100644 --- a/lisp/interface/matrix.lisp +++ b/lisp/interface/matrix.lisp @@ -120,3 +120,7 @@ MATRIX.")) (values start-row start-column end-row end-column) (error "The matrix range (~D:~D,~D:~D) is invalid." start-row start-column end-row end-column)))) + +(defgeneric mat-equal (matrix1 matrix2) + (:documentation + "Return if the two arrays elments are the same")) diff --git a/lisp/interface/vector.lisp b/lisp/interface/vector.lisp index 67a5d31..51cb9bb 100644 --- a/lisp/interface/vector.lisp +++ b/lisp/interface/vector.lisp @@ -89,6 +89,10 @@ "Destructively modifies the result vector with the result of applying the function to each element of the vectors.")) +(defgeneric reduce-vector (function vector initial-value) + (:documentation + "Reduce a vector into a scalar using a function")) + (defmacro dovector ((element vector &optional result) &body body) "Iterate over vector returning result." (let ((pos (gensym "POS-")) @@ -108,3 +112,15 @@ applying the function to each element of the vectors.")) (defgeneric napply-rotation (vector1 vector2 cc ss) (:documentation "Return the plane rotations of vector1 and vector2 by cc and ss.")) + +(defgeneric vec-equal (vector1 vector2) + (:documentation + "Test if the contents of vector1 and vector2 are the same")) + +(defgeneric elem-divide (vector1 vector2) + (:documentation + "Return the element by element division of two vectors")) + +(defgeneric elem-multiply (vector1 vector2) + (:documentation + "Return the element by element multiplation of two vectors")) diff --git a/lisp/kernel/binary-operations.lisp b/lisp/kernel/binary-operations.lisp index 5cfa579..f4948e4 100644 --- a/lisp/kernel/binary-operations.lisp +++ b/lisp/kernel/binary-operations.lisp @@ -75,6 +75,16 @@ the operation.")) "Return the scaled operation." (lambda (n1 n2) (- (* scalar1 n1) (* scalar2 n2)))) +(defmethod binary-op + ((op (eql #'/)) ) + "Return the scaled operation." + (lambda (n1 n2) (/ n1 n2))) + +(defmethod binary-op + ((op (eql #'*)) ) + "Return the scaled operation." + (lambda (n1 n2) (* n1 n2))) + ;;; Binary vector operations (defun %vector<-vector1-op-vector2 (operation vector1 vector2) @@ -124,6 +134,18 @@ addition." (scaled-binary-op #'- scalar1 scalar2) vector1 vector2)) +(defun element-divide-vector (vector1 vector2) + "Vector pointwise division." + (%vector<-vector1-op-vector2 + (binary-op #'/) + vector1 vector2)) + +(defun element-multiply-vector (vector1 vector2) + "Vector pointwise multiplication." + (%vector<-vector1-op-vector2 + (binary-op #'*) + vector1 vector2)) + (defun nsubtract-vector (vector1 vector2 scalar1 scalar2) "Destructive vector binary subtraction." (%vector1<-vector1-op-vector2 @@ -137,7 +159,7 @@ addition." and element2 across vector2 sum (* element1 element2) into result finally - (return (if scalar (* scalar result) result)))) + (return (if scalar (* scalar result) result)))) ;;; Binary array/vector operations diff --git a/lisp/kernel/linear-algebra-kernel.lisp b/lisp/kernel/linear-algebra-kernel.lisp index 54b1bd0..c3b082c 100644 --- a/lisp/kernel/linear-algebra-kernel.lisp +++ b/lisp/kernel/linear-algebra-kernel.lisp @@ -52,6 +52,8 @@ :subtract-vector :nsubtract-vector :add-array :nadd-array :subtract-array :nsubtract-array + :element-multiply-vector + :element-divide-vector :inner-product-vector :product-vector-array :product-array-vector diff --git a/lisp/linear-algebra.lisp b/lisp/linear-algebra.lisp index 4062d58..138fd24 100644 --- a/lisp/linear-algebra.lisp +++ b/lisp/linear-algebra.lisp @@ -51,7 +51,12 @@ :replace-vector :map-vector :map-into-vector + :reduce-vector :dovector + :vec-equal + :elem-divide + :elem-multiply + :distance :apply-rotation :napply-rotation) ;; Matrix interface (:export :matrix-object From e027f627b7e10113aedf5db9f378a8882d8fe9a2 Mon Sep 17 00:00:00 2001 From: Brian Eberman Date: Sun, 4 Sep 2022 14:11:27 -0400 Subject: [PATCH 04/12] added a distance function and a print function --- lisp/data-vector.lisp | 8 ++++++++ lisp/interface/vector.lisp | 5 +++++ lisp/linear-algebra.lisp | 3 ++- 3 files changed, 15 insertions(+), 1 deletion(-) diff --git a/lisp/data-vector.lisp b/lisp/data-vector.lisp index 97245c8..ed6678f 100644 --- a/lisp/data-vector.lisp +++ b/lisp/data-vector.lisp @@ -36,6 +36,10 @@ (:documentation "A data vector.")) +(defmethod print-object ((vector data-vector) stream) + (print-unreadable-object (vector stream :type t :identity t) + (format stream "~a" (contents vector)))) + (defclass row-vector (data-vector) () (:documentation @@ -494,6 +498,10 @@ applying the function to each element of the vectors." (contents vector1) (contents vector2)) )) +(defmethod vec-every ((vector data-vector) predicate) + (every predicate (contents vector))) + + (defmethod distance :before ((vector1 data-vector) (vector2 data-vector) &optional (measure 2)) "Verify that the dimensions are equal." diff --git a/lisp/interface/vector.lisp b/lisp/interface/vector.lisp index 51cb9bb..6ed6478 100644 --- a/lisp/interface/vector.lisp +++ b/lisp/interface/vector.lisp @@ -124,3 +124,8 @@ applying the function to each element of the vectors.")) (defgeneric elem-multiply (vector1 vector2) (:documentation "Return the element by element multiplation of two vectors")) + +(defgeneric vec-every (vector predicate) + (:documentation + "Return true if every element in the vector satisfies + predicate. Return nil on the first element that does not")) diff --git a/lisp/linear-algebra.lisp b/lisp/linear-algebra.lisp index 138fd24..58e89be 100644 --- a/lisp/linear-algebra.lisp +++ b/lisp/linear-algebra.lisp @@ -48,7 +48,7 @@ :vref :copy-vector :subvector - :replace-vector + :replace-vectorx :map-vector :map-into-vector :reduce-vector @@ -56,6 +56,7 @@ :vec-equal :elem-divide :elem-multiply + :vec-everyx :distance :apply-rotation :napply-rotation) ;; Matrix interface From 3c8b128d2043cedcbc5bd35ce4b95efa713a5ea1 Mon Sep 17 00:00:00 2001 From: Brian Eberman Date: Sat, 19 Nov 2022 09:33:09 -0500 Subject: [PATCH 05/12] Added a min function on vector which finds the smallest absolute distance between two vectorl a make-identity-matrix to create an identity matrix; a norm of vectors --- lisp/data-vector.lisp | 4 ++++ lisp/interface/fundamental-ops.lisp | 4 ++++ lisp/interface/identity-matrix.lisp | 5 +++++ lisp/kernel/unary-operations.lisp | 2 ++ lisp/linear-algebra.lisp | 8 ++++++-- 5 files changed, 21 insertions(+), 2 deletions(-) diff --git a/lisp/data-vector.lisp b/lisp/data-vector.lisp index ed6678f..34cea8c 100644 --- a/lisp/data-vector.lisp +++ b/lisp/data-vector.lisp @@ -267,6 +267,10 @@ applying the function to each element of the vectors." "Return the p-norm of the vector." (norm-vector (contents vector) measure)) +(defmethod min-vector ((vector data-vector)) + (loop for element across (contents vector) minimize (abs element))) + + (defmethod transpose ((vector column-vector)) "Return a row vector." (make-instance diff --git a/lisp/interface/fundamental-ops.lisp b/lisp/interface/fundamental-ops.lisp index f438f71..69e4b8a 100644 --- a/lisp/interface/fundamental-ops.lisp +++ b/lisp/interface/fundamental-ops.lisp @@ -32,6 +32,10 @@ (:documentation "Return the norm according to measure.")) +(defgeneric min-vector (vector) + (:documentation + "Return the minimum absolute value")) + (defgeneric distance (vector1 vector2 &optional measure) (:documentation "Return the distance between vector1 and vector2 using the measure")) diff --git a/lisp/interface/identity-matrix.lisp b/lisp/interface/identity-matrix.lisp index 83ba84e..0c2d32e 100644 --- a/lisp/interface/identity-matrix.lisp +++ b/lisp/interface/identity-matrix.lisp @@ -63,6 +63,11 @@ :initial-contents (list (coerce 0 element-type) (coerce 1 element-type))))))) +(defmethod make-identity-matrix (&key dimension element-type ) + (make-instance 'identity-matrix + :dimensions (list dimension dimension) :element-type element-type)) + + (defmethod matrix-in-bounds-p ((matrix identity-matrix) (row integer) (column integer)) "Return true if row and column do not exceed the dimensions of matrix." diff --git a/lisp/kernel/unary-operations.lisp b/lisp/kernel/unary-operations.lisp index 3421c20..2eedc29 100644 --- a/lisp/kernel/unary-operations.lisp +++ b/lisp/kernel/unary-operations.lisp @@ -193,6 +193,7 @@ array column." (dotimes (index (length vector) result) (setf (aref result index) (abs (aref vector index)))))) + (defmethod norm-vector ((data vector) (measure (eql 1))) "Return the Taxicab norm of the list." (loop for element across data sum (abs element))) @@ -213,6 +214,7 @@ array column." "Return the infinity, or maximum, norm of vector." (loop for element across data maximize (abs element))) + (defmethod norm-array ((data array) (measure (eql 1))) "Return the 1 norm of the array." (let ((m-rows (array-dimension data 0)) diff --git a/lisp/linear-algebra.lisp b/lisp/linear-algebra.lisp index 58e89be..1b64e99 100644 --- a/lisp/linear-algebra.lisp +++ b/lisp/linear-algebra.lisp @@ -31,6 +31,7 @@ (:use :floating-point :linear-algebra-kernel) ;; Fundamental operations (:export :norm + :min-vector :transpose :ntranspose :permute :scale :nscale @@ -76,8 +77,11 @@ :replace-matrix :matrix-validated-range) ;; Identity matrix - (:export :identity-matrix - :identity-matrix-p) + (:export + :identity-matrix + :identity-matrix-p + :make-identity-matrix + ) ;; Permutation matrix (:export :permutation-matrix :permutation-matrix-p) From e610d5af56f74ac291f98dd93d7e96348001ff1d Mon Sep 17 00:00:00 2001 From: Brian Eberman Date: Sun, 20 Nov 2022 09:27:48 -0500 Subject: [PATCH 06/12] Fixed a few tests. Added the function add-diagonal --- lisp/dense-matrix.lisp | 28 ++++++++++++++++++++++ lisp/interface/fundamental-ops.lisp | 4 ++++ lisp/interface/identity-matrix.lisp | 5 ---- lisp/linear-algebra.lisp | 7 +++--- test/data-vector.lisp | 15 ++++++------ test/dense-matrix.lisp | 37 +++++++++++++++++++++++++++++ test/linear-algebra-test.lisp | 1 + 7 files changed, 82 insertions(+), 15 deletions(-) diff --git a/lisp/dense-matrix.lisp b/lisp/dense-matrix.lisp index 4d3d1c8..35b234c 100644 --- a/lisp/dense-matrix.lisp +++ b/lisp/dense-matrix.lisp @@ -495,3 +495,31 @@ matrix with a column vector." :contents (gauss-invert (contents matrix))) (error "The number of rows does not equal columns."))) + +(defmethod add-diagonal ((value number) (matrix dense-matrix)) + (let ((result (copy-matrix matrix)) + (num (min (matrix-dimensions matrix))) + ) + (dotimes (index num result) + (setf + (aref result index index) + (+ (aref matrix index index) value))) + result)) + +(defmethod add-diagonal :before + ((vec data-vector) (matrix dense-matrix)) + (unless (compatible-dimensions-p :solve matrix vector) + (error "Matrix~A is incompatible with column vector(~D)." + (matrix-dimensions matrix) (vector-length vector)))) + +(defmethod add-diagonal ((vec data-vector) (matrix dense-matrix)) + (let ((result (copy-matrix matrix)) + (num (min (matrix-dimensions matrix))) + ) + (dotimes (index num result) + (setf + (aref result index index) + (+ (aref matrix index index) + (vref vec index)))) + result)) + diff --git a/lisp/interface/fundamental-ops.lisp b/lisp/interface/fundamental-ops.lisp index 69e4b8a..1716c3f 100644 --- a/lisp/interface/fundamental-ops.lisp +++ b/lisp/interface/fundamental-ops.lisp @@ -101,5 +101,9 @@ (:documentation "Return the invert of the matrix with in-place decomposition.")) +(defgeneric add-diagonal (scalar-or-vector matrix) + (:documentation + "Add a scalar or vector to every diagonal of the matrix")) + diff --git a/lisp/interface/identity-matrix.lisp b/lisp/interface/identity-matrix.lisp index 0c2d32e..83ba84e 100644 --- a/lisp/interface/identity-matrix.lisp +++ b/lisp/interface/identity-matrix.lisp @@ -63,11 +63,6 @@ :initial-contents (list (coerce 0 element-type) (coerce 1 element-type))))))) -(defmethod make-identity-matrix (&key dimension element-type ) - (make-instance 'identity-matrix - :dimensions (list dimension dimension) :element-type element-type)) - - (defmethod matrix-in-bounds-p ((matrix identity-matrix) (row integer) (column integer)) "Return true if row and column do not exceed the dimensions of matrix." diff --git a/lisp/linear-algebra.lisp b/lisp/linear-algebra.lisp index 1b64e99..3f4402a 100644 --- a/lisp/linear-algebra.lisp +++ b/lisp/linear-algebra.lisp @@ -49,7 +49,7 @@ :vref :copy-vector :subvector - :replace-vectorx + :replace-vector :map-vector :map-into-vector :reduce-vector @@ -75,12 +75,13 @@ :mat-equal :submatrix :replace-matrix - :matrix-validated-range) + :matrix-validated-range + :add-diagonal + ) ;; Identity matrix (:export :identity-matrix :identity-matrix-p - :make-identity-matrix ) ;; Permutation matrix (:export :permutation-matrix diff --git a/test/data-vector.lisp b/test/data-vector.lisp index 9da481f..d548025 100644 --- a/test/data-vector.lisp +++ b/test/data-vector.lisp @@ -1207,18 +1207,18 @@ (linear-algebra:product (linear-algebra:row-vector 1 2 3) (linear-algebra:column-vector 1 2 3 4))) - (assert-error - 'error + (assert-equal + 14 (linear-algebra:product (linear-algebra:column-vector 1 2 3) (linear-algebra:column-vector 1 2 3))) - (assert-error - 'error + (assert-equal + 14 (linear-algebra:product (linear-algebra:row-vector 1 2 3) (linear-algebra:row-vector 1 2 3))) - (assert-error - 'error + (assert-equal + 14 (linear-algebra:product (linear-algebra:column-vector 1 2 3) (linear-algebra:row-vector 1 2 3)))) @@ -1231,5 +1231,6 @@ 14 (linear-algebra-kernel:inner-product-vector (linear-algebra:row-vector 1 2 3) - (linear-algebra:row-vector 1 2 3))) + (linear-algebra:row-vector 1 2 3) + nil)) ) diff --git a/test/dense-matrix.lisp b/test/dense-matrix.lisp index 5307109..f3b77df 100644 --- a/test/dense-matrix.lisp +++ b/test/dense-matrix.lisp @@ -1363,3 +1363,40 @@ ( 0.031 0.015815994 29.992375 0.33191067) ( 0.041 0.020756614 0.014001981 39.98469)) matrix))) + +(define-test add-diagonal-dense-matrix + (:tag :dense-matrix :add-diagonal) + (let ((*epsilon* (* 64 single-float-epsilon)) + (scalar 2.0) + (vector1 (linear-algebra:column-vector 1.0 2.0)) + (vector2 (linear-algebra:column-vector 1.0 2.0 3.0)) + (matrix1 + (linear-algebra:make-matrix + 2 2 :initial-contents '((1.1 1.2) (2.1 2.2)))) + (result1 + (linear-algebra:make-matrix + 2 2 :initial-contents + '((3.1 1.2) (2.1 4.2)))) + (result2 + (linear-algebra:make-matrix + 2 2 :initial-contents + '((2.1 1.2) (2.1 4.2)))) + ) + + ;; 2x2 + (assert-true + (linear-algebra:mat-equal + result1 + (linear-algebra:add-diagonal scalar matrix1))) + + (assert-true + (linear-algebra:mat-equal + result2 + (linear-algebra:add-diagonal vector1 matrix1))) + + (assert-error + 'error + (linear-algebra:add-diagonal vector2 matrix1)) + )) + + diff --git a/test/linear-algebra-test.lisp b/test/linear-algebra-test.lisp index 4b34341..ee62e2f 100644 --- a/test/linear-algebra-test.lisp +++ b/test/linear-algebra-test.lisp @@ -7,6 +7,7 @@ (in-package :linear-algebra-test) + ;;; Convenience functions (defun random-interior-index (size) From 5944cf382f49fca61e5435073f7c87f4ec64f27a Mon Sep 17 00:00:00 2001 From: Brian Eberman Date: Sun, 20 Nov 2022 09:39:23 -0500 Subject: [PATCH 07/12] adding a todo and claritying the README --- README | 44 ++++++++++++++++++++++++++++++++++++++++++++ todo.org | 11 +++++++++++ 2 files changed, 55 insertions(+) create mode 100644 todo.org diff --git a/README b/README index c883e80..da90d01 100644 --- a/README +++ b/README @@ -1,6 +1,50 @@ +Linear Algebra in Common Lisp + +Copyright (c) 2011-2014, Odonata Research LLC + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. + +Copyright (c) 2021-2023, Ten Factor Growth, LLC +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. + + linear-algebra linear-algebra is a Common Lisp library of numeric linear algebra routines. It approximates the functionality present in the Basic Linear Algebra Subroutines (BLAS) Fortran library. The objective of this project is to generate linear algebra routines in native ANSI Common Lisp. + +Use within quicklisp with (ql:quickload :linear-algebra) diff --git a/todo.org b/todo.org new file mode 100644 index 0000000..f4dcfd6 --- /dev/null +++ b/todo.org @@ -0,0 +1,11 @@ +* TODO convert the tests to fiveam +* TODO move the types into a types subdirectory +** create a types +** refactor dense-matrix to remove the operations and put them in fundamental-ops +* TODO get identity matrix to work +** put it in types +** create a diagonal matrix type and make identity a subclass +* TODO add banded matrix type +* TODO add SVD decomposition for Hermetian +* TODO add gram-schmidt QR +* TODO add sequential gram-schmidt as a square root op From e25abbe7b01c5fae367cc58135fc36d19f1282ff Mon Sep 17 00:00:00 2001 From: Brian Eberman Date: Sun, 20 Nov 2022 11:06:23 -0500 Subject: [PATCH 08/12] fixed the add-diagonal --- lisp/dense-matrix.lisp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/lisp/dense-matrix.lisp b/lisp/dense-matrix.lisp index 35b234c..32a0950 100644 --- a/lisp/dense-matrix.lisp +++ b/lisp/dense-matrix.lisp @@ -498,12 +498,12 @@ matrix with a column vector." (defmethod add-diagonal ((value number) (matrix dense-matrix)) (let ((result (copy-matrix matrix)) - (num (min (matrix-dimensions matrix))) + (num (apply #'min (matrix-dimensions matrix))) ) (dotimes (index num result) (setf - (aref result index index) - (+ (aref matrix index index) value))) + (mref result index index) + (+ (mref matrix index index) value))) result)) (defmethod add-diagonal :before @@ -518,8 +518,8 @@ matrix with a column vector." ) (dotimes (index num result) (setf - (aref result index index) - (+ (aref matrix index index) + (mref result index index) + (+ (mref matrix index index) (vref vec index)))) result)) From b65654cdb96dfbbace709038e1a6a2380bbcaa0d Mon Sep 17 00:00:00 2001 From: Brian Eberman Date: Tue, 17 Jan 2023 09:53:05 -0500 Subject: [PATCH 09/12] squash! fixed the add-diagonal --- lisp/data-vector.lisp | 5 +++++ lisp/dense-matrix.lisp | 3 +++ lisp/interface/vector.lisp | 4 ++++ lisp/kernel/binary-operations.lisp | 9 +++++++++ lisp/kernel/linear-algebra-kernel.lisp | 1 + lisp/linear-algebra.lisp | 3 ++- 6 files changed, 24 insertions(+), 1 deletion(-) diff --git a/lisp/data-vector.lisp b/lisp/data-vector.lisp index 34cea8c..c3000f5 100644 --- a/lisp/data-vector.lisp +++ b/lisp/data-vector.lisp @@ -502,6 +502,11 @@ applying the function to each element of the vectors." (contents vector1) (contents vector2)) )) +(defmethod elem-greater + ((vector1 data-vector) (vector2 data-vector)) + (element-greater-vector + (contents vector1) (contents vector2))) + (defmethod vec-every ((vector data-vector) predicate) (every predicate (contents vector))) diff --git a/lisp/dense-matrix.lisp b/lisp/dense-matrix.lisp index 32a0950..5354066 100644 --- a/lisp/dense-matrix.lisp +++ b/lisp/dense-matrix.lisp @@ -34,6 +34,9 @@ (:documentation "Dense matrix object.")) +(defmethod print-object ((matrix dense-matrix) stream) + (print-unreadable-object (matrix stream :type t :identity t) + (format stream "~a" (contents matrix)))) (defmethod initialize-instance :after ((self dense-matrix) &rest initargs diff --git a/lisp/interface/vector.lisp b/lisp/interface/vector.lisp index 6ed6478..c956d39 100644 --- a/lisp/interface/vector.lisp +++ b/lisp/interface/vector.lisp @@ -125,6 +125,10 @@ applying the function to each element of the vectors.")) (:documentation "Return the element by element multiplation of two vectors")) +(defgeneric elem-greater (vector1 vector2) + (:documentation + "Return if vector1 is greater than vector 2 on every element")) + (defgeneric vec-every (vector predicate) (:documentation "Return true if every element in the vector satisfies diff --git a/lisp/kernel/binary-operations.lisp b/lisp/kernel/binary-operations.lisp index f4948e4..2dec6d4 100644 --- a/lisp/kernel/binary-operations.lisp +++ b/lisp/kernel/binary-operations.lisp @@ -146,6 +146,15 @@ addition." (binary-op #'*) vector1 vector2)) +(defun element-greater-vector (vector1 vector2) + (loop + for element1 across vector1 + and element2 across vector2 + when (<= element1 element2) + return nil + finally (return t))) + + (defun nsubtract-vector (vector1 vector2 scalar1 scalar2) "Destructive vector binary subtraction." (%vector1<-vector1-op-vector2 diff --git a/lisp/kernel/linear-algebra-kernel.lisp b/lisp/kernel/linear-algebra-kernel.lisp index c3b082c..5000233 100644 --- a/lisp/kernel/linear-algebra-kernel.lisp +++ b/lisp/kernel/linear-algebra-kernel.lisp @@ -54,6 +54,7 @@ :subtract-array :nsubtract-array :element-multiply-vector :element-divide-vector + :element-greater-vector :inner-product-vector :product-vector-array :product-array-vector diff --git a/lisp/linear-algebra.lisp b/lisp/linear-algebra.lisp index 3f4402a..3e00546 100644 --- a/lisp/linear-algebra.lisp +++ b/lisp/linear-algebra.lisp @@ -57,7 +57,8 @@ :vec-equal :elem-divide :elem-multiply - :vec-everyx + :vec-every + :elem-greater :distance :apply-rotation :napply-rotation) ;; Matrix interface From 3b7a228d6b80d0118211b802884609b9c161dcae Mon Sep 17 00:00:00 2001 From: Brian Eberman Date: Tue, 31 Jan 2023 20:38:10 -0500 Subject: [PATCH 10/12] Added a distance function --- lisp/data-vector.lisp | 2 ++ lisp/interface/fundamental-ops.lisp | 2 -- lisp/interface/vector.lisp | 1 + 3 files changed, 3 insertions(+), 2 deletions(-) diff --git a/lisp/data-vector.lisp b/lisp/data-vector.lisp index c3000f5..0900fe4 100644 --- a/lisp/data-vector.lisp +++ b/lisp/data-vector.lisp @@ -521,3 +521,5 @@ applying the function to each element of the vectors." (defmethod distance ((vector1 data-vector) (vector2 data-vector) &optional (measure 2)) (norm (subtract vector1 vector2) measure)) + + diff --git a/lisp/interface/fundamental-ops.lisp b/lisp/interface/fundamental-ops.lisp index 1716c3f..44fabe6 100644 --- a/lisp/interface/fundamental-ops.lisp +++ b/lisp/interface/fundamental-ops.lisp @@ -105,5 +105,3 @@ (:documentation "Add a scalar or vector to every diagonal of the matrix")) - - diff --git a/lisp/interface/vector.lisp b/lisp/interface/vector.lisp index c956d39..6f83bd6 100644 --- a/lisp/interface/vector.lisp +++ b/lisp/interface/vector.lisp @@ -133,3 +133,4 @@ applying the function to each element of the vectors.")) (:documentation "Return true if every element in the vector satisfies predicate. Return nil on the first element that does not")) + From e35cb74dcf36549afc34c85cfbb3827ac49877e5 Mon Sep 17 00:00:00 2001 From: Brian Eberman Date: Tue, 31 Jan 2023 21:32:48 -0500 Subject: [PATCH 11/12] added matrix-trace and outer-product --- lisp/data-vector.lisp | 24 +++++++++++++++++++++++- lisp/dense-matrix.lisp | 13 +++++++++++++ lisp/interface/matrix.lisp | 4 ++++ lisp/interface/vector.lisp | 3 +++ lisp/linear-algebra.lisp | 7 +++++-- 5 files changed, 48 insertions(+), 3 deletions(-) diff --git a/lisp/data-vector.lisp b/lisp/data-vector.lisp index 0900fe4..f9756f2 100644 --- a/lisp/data-vector.lisp +++ b/lisp/data-vector.lisp @@ -522,4 +522,26 @@ applying the function to each element of the vectors." ((vector1 data-vector) (vector2 data-vector) &optional (measure 2)) (norm (subtract vector1 vector2) measure)) - +(defmethod outer-product ((vector1 data-vector) + (vector2 data-vector)) + (let* ((len1 (vector-length vector1)) + (len2 (vector-length vector2)) + (mat + (if (equal vector1 vector2) + (make-instance 'square-matrix + :dimensions (list len1 len2) + :initial-element 0 + :element-type (vector-element-type vector1)) + (make-instance 'dense-matrix + :dimensions (list len1 len2) + :initial-element 0 + :element-type (vector-element-type vector1)))) + ) + (dotimes (i len1) + (dotimes (j len2) + (setf + (mref mat i j) + (* (vref vector1 i) (vref vector2 j))) + )) + mat + )) diff --git a/lisp/dense-matrix.lisp b/lisp/dense-matrix.lisp index 5354066..975d9f5 100644 --- a/lisp/dense-matrix.lisp +++ b/lisp/dense-matrix.lisp @@ -526,3 +526,16 @@ matrix with a column vector." (vref vec index)))) result)) +(defmethod matrix-trace ((matrix dense-matrix)) + (let* ((num (apply #'min (matrix-dimensions matrix))) + (vec (make-instance 'data-vector :size num + :initial-element 0 + :element-type (matrix-element-type matrix))) + ) + (dotimes (index num ) + (setf + (vref vec index) + (mref matrix index index)) + ) + vec + )) diff --git a/lisp/interface/matrix.lisp b/lisp/interface/matrix.lisp index f01e9e7..f4cf9a1 100644 --- a/lisp/interface/matrix.lisp +++ b/lisp/interface/matrix.lisp @@ -99,6 +99,10 @@ MATRIX.")) (:documentation "Set the submatrix of the matrix.")) +(defgeneric matrix-trace (matrix) + (:documentation + "Return a column vector that is the trace of matrix")) + (defgeneric replace-matrix (matrix1 matrix2 &key start-row1 end-row1 diff --git a/lisp/interface/vector.lisp b/lisp/interface/vector.lisp index 6f83bd6..e802e29 100644 --- a/lisp/interface/vector.lisp +++ b/lisp/interface/vector.lisp @@ -134,3 +134,6 @@ applying the function to each element of the vectors.")) "Return true if every element in the vector satisfies predicate. Return nil on the first element that does not")) +(defgeneric outer-product (vec1 vec2) + (:documentation + "Outer product of vec1 and vec2 returns an nxm matrix")) diff --git a/lisp/linear-algebra.lisp b/lisp/linear-algebra.lisp index 3e00546..2a375ab 100644 --- a/lisp/linear-algebra.lisp +++ b/lisp/linear-algebra.lisp @@ -58,9 +58,11 @@ :elem-divide :elem-multiply :vec-every + :outer-product :elem-greater :distance - :apply-rotation :napply-rotation) + :apply-rotation :napply-rotation + ) ;; Matrix interface (:export :matrix-object :initialize-matrix @@ -73,7 +75,8 @@ :matrix-column-dimension :mref :copy-matrix - :mat-equal + :mat-equal + :trace :submatrix :replace-matrix :matrix-validated-range From 279a2f40aec562a40aed3c4babd782c0f885c40a Mon Sep 17 00:00:00 2001 From: Brian Eberman Date: Wed, 1 Feb 2023 08:47:10 -0500 Subject: [PATCH 12/12] fixed trace to call it diagonal since the function returns the diagonal. Enhanced outer-product so that it determines the appropriate matrix type based on the two vectors. --- lisp/data-vector.lisp | 18 ++++++++---------- lisp/dense-matrix.lisp | 2 +- lisp/interface/matrix.lisp | 4 ++++ lisp/linear-algebra.lisp | 4 ++-- 4 files changed, 15 insertions(+), 13 deletions(-) diff --git a/lisp/data-vector.lisp b/lisp/data-vector.lisp index f9756f2..3d39391 100644 --- a/lisp/data-vector.lisp +++ b/lisp/data-vector.lisp @@ -526,16 +526,14 @@ applying the function to each element of the vectors." (vector2 data-vector)) (let* ((len1 (vector-length vector1)) (len2 (vector-length vector2)) - (mat - (if (equal vector1 vector2) - (make-instance 'square-matrix - :dimensions (list len1 len2) - :initial-element 0 - :element-type (vector-element-type vector1)) - (make-instance 'dense-matrix - :dimensions (list len1 len2) - :initial-element 0 - :element-type (vector-element-type vector1)))) + (mat-type (cond + ((equal vector1 vector2) 'symmetric-matrix) + ((= len1 len2) 'square-matrix) + (t 'dense-matrix))) + (mat (make-instance mat-type + :dimensions (list len1 len2) + :initial-element 0 + :element-type (vector-element-type vector1))) ) (dotimes (i len1) (dotimes (j len2) diff --git a/lisp/dense-matrix.lisp b/lisp/dense-matrix.lisp index 975d9f5..73e088d 100644 --- a/lisp/dense-matrix.lisp +++ b/lisp/dense-matrix.lisp @@ -526,7 +526,7 @@ matrix with a column vector." (vref vec index)))) result)) -(defmethod matrix-trace ((matrix dense-matrix)) +(defmethod matrix-diagonal ((matrix dense-matrix)) (let* ((num (apply #'min (matrix-dimensions matrix))) (vec (make-instance 'data-vector :size num :initial-element 0 diff --git a/lisp/interface/matrix.lisp b/lisp/interface/matrix.lisp index f4cf9a1..1eb72d3 100644 --- a/lisp/interface/matrix.lisp +++ b/lisp/interface/matrix.lisp @@ -128,3 +128,7 @@ MATRIX.")) (defgeneric mat-equal (matrix1 matrix2) (:documentation "Return if the two arrays elments are the same")) + +(defgeneric matrix-diagonal (matrix) + (:documentation + "Return the diagonal of the matrix as a vector")) diff --git a/lisp/linear-algebra.lisp b/lisp/linear-algebra.lisp index 2a375ab..b31f16b 100644 --- a/lisp/linear-algebra.lisp +++ b/lisp/linear-algebra.lisp @@ -61,7 +61,7 @@ :outer-product :elem-greater :distance - :apply-rotation :napply-rotation + :apply-rotation :napply-rotation ) ;; Matrix interface (:export :matrix-object @@ -76,7 +76,7 @@ :mref :copy-matrix :mat-equal - :trace + :matrix-diagonal :submatrix :replace-matrix :matrix-validated-range