From 1e5bd820f87349cd7cfb0411edb4182ca4d89119 Mon Sep 17 00:00:00 2001 From: Cat Stevens Date: Sat, 29 Feb 2020 21:46:38 -0500 Subject: [PATCH] math.matrices: fix/rename mnorm, update all norms closes #2244 - `mnorm` has been renamed to `normalize-matrix` to reflect what it actually does, which is normalize a matrix, not find a norm of a matrix. - `mnorm` is no longer a word defined here. - bugfix: previously, `normalize-matrix` found the supremum of a matrix (`mmax`), before taking the supremum's absolute value (`abs`) and dividing the matrix by it (`m/n`). for matrices containing only negative values and 0, the supremum is 0, and a `div-by-zero` error was thrown. `normalize-matrix` has been fixed to first `abs` all the matrix elements, and then find the supremum and divide, it also receieved a zero-matrix? guard for optimization and preventing `div-by-zero`. - new alias: `hilbert-schmidt-norm` for `frobenius-norm`, to go along with `math.matrices.extras.` and improve searchability by physicists. - new word: `matrix-p-norm`, written as an analogue of `math.vectors.p-norm`. - new word: `matrix-p-q-norm`, which generalizes entrywise matrix norm over the L^p,q vector space. - new word: `matrix-p-norm-entrywise`: `matrix-p-norm`'s fallback for p =/= 1, 2, inf; analogue of `math.vectors.p-norm-default`. - all norm words have gotten new docs, `zero-matrix?` guards as an optimisation, and most have gotten new tests. --- basis/math/matrices/matrices-docs.factor | 136 ++++++++++++++---- basis/math/matrices/matrices-tests.factor | 72 ++++------ basis/math/matrices/matrices.factor | 52 ++++++- extra/math/matrices/extras/extras-docs.factor | 3 - 4 files changed, 186 insertions(+), 77 deletions(-) diff --git a/basis/math/matrices/matrices-docs.factor b/basis/math/matrices/matrices-docs.factor index 9bf0da32fa..cd06a673a6 100644 --- a/basis/math/matrices/matrices-docs.factor +++ b/basis/math/matrices/matrices-docs.factor @@ -1,4 +1,4 @@ -! Copyright (C) 2005, 2010, 2018 Slava Pestov, Joe Groff, and Cat Stevens. +! Copyright (C) 2005, 2010, 2018, 2020 Slava Pestov, Joe Groff, and Cat Stevens. USING: accessors arrays assocs generic.single formatting locals help.markup help.markup.private help.syntax io kernel math math.functions math.order math.ratios math.vectors opengl.gl prettyprint sequences sequences.generalizations urls ; @@ -136,7 +136,15 @@ $nl } { $subsections main-diagonal anti-diagonal +} +"The following matrix norms are provided in the 𝑙ₚ vector space; these words are equivalent to ∥･∥ₚ for " { $snippet "p = 1, 2, ∞, n" } ", respectively:" +{ $subsections + m-1norm + m-infinity-norm + frobenius-norm + matrix-p-norm-entrywise + matrix-p-norm } ; ! PREDICATE CLASSES @@ -524,6 +532,17 @@ HELP: cols } } ; +HELP: >square-matrix +{ $values { "m" matrix } { "subset" square-matrix } } +{ $description "Find only the " { $link2 square-matrix "square" } " subset of the input matrix." } +{ $examples + { $example + "USING: math.matrices prettyprint ;" + "{ { 0 2 4 6 } { 1 3 5 7 } } >square-matrix ." + "{ { 0 2 } { 1 3 } }" + } +} ; + HELP: matrix-map { $values { "matrix" matrix } { "quot" { $quotation ( ... elt -- ... elt' ) } } { "matrix'" matrix } } { $description "Apply the quotation to every element of the matrix." } @@ -938,7 +957,97 @@ HELP: mmax } } ; -HELP: mnorm +HELP: m-1norm +{ $values { "m" matrix } { "n" number } } +{ $description "Find the size of a matrix in 𝑙₁ (" { $snippet "L^₁" } ") vector space, usually written ∥･∥₁." +$nl "This is the matrix norm when " { $snippet "p=1" } ", and is the overall maximum of the sums of the columns." } +{ $notelist { $equiv-word-note "transpose" m-infinity-norm } } +{ $examples + { $example + "USING: math.matrices prettyprint ;" + "{ { 2 -2 1 } { 1 3 -1 } { 2 -4 2 } } m-1norm ." + "9" + } +} ; + +HELP: m-infinity-norm +{ $values { "m" matrix } { "n" number } } +{ $description "Find the size of a matrix, in 𝑙∞ (" { $snippet "L^∞" } ") vector space, usually written ∥･∥∞." +$nl "This is the matrix norm when " { $snippet "p=∞" } ", and is the overall maximum of the sums of the rows." } +{ $notelist { $equiv-word-note "transpose" m-1norm } } +{ $examples + { $example + "USING: math.matrices prettyprint ;" + "{ { 2 -2 1 } { 1 3 -1 } { 2 -4 2 } } m-infinity-norm ." + "8" + } +} ; + +HELP: frobenius-norm +{ $values { "m" matrix } { "n" number } } +{ $description "Find the size of a matrix in 𝑙₂ (" { $snippet "L^2" } ") vector space, usually written ∥･∥₂ₚ." +$nl "This is the matrix norm when " { $snippet "p=2" } ", and is the square root of the sums of the squares of all the elements of the matrix." } +{ $notelist + { "This norm is sometimes called the Hilbert-Schmidt norm." } + { "Because " $snippet { "p=2" } ", this word could be named " { $snippet "m-2norm" } "." } +} +{ $examples + { $example + "USING: math.matrices prettyprint ;" + "{ { 1 1 } { 1 1 } } frobenius-norm ." + "2.0" + } +} ; + +{ m-1norm m-infinity-norm frobenius-norm } related-words + +HELP: matrix-p-q-norm +{ $values { "m" matrix } { "p" "positive real number" } { "q" "positive real number" } { "n" "non-negative real number" } } +{ $description "Find the size of a matrix in " { $snippet "L^p,q" } " vector space." +$nl "This is the matrix norm for any " { $snippet "p, q ≥ 1" } ". It is still an entry-wise norm, like " { $link matrix-p-norm-entrywise } "." } +{ $examples + "Equivalent to " { $link frobenius-norm } " for " { $snippet "p = q = 2 " } ":" + { $example + "USING: math.matrices prettyprint ;" + "{ { 1 1 } { 1 1 } } 2 2 matrix-p-q-norm ." + "2.0" + } +} ; + +HELP: matrix-p-norm-entrywise +{ $values { "m" matrix } { "p" "positive real number" } { "n" "non-negative real number" } } +{ $description "Find the entry-wise norm of a matrix, in 𝑙ₚ (" { $snippet "L^p" } ") vector space." } +{ $notes "This word is distinct from a Schatten p-norm, as well as any of " { $links m-1norm frobenius-norm m-infinity-norm } "." } +{ $examples + { $example + "USING: math.matrices prettyprint ;" + "4 4 1 2 matrix-p-norm-entrywise ." + "4.0" + } +} ; + +HELP: matrix-p-norm +{ $values { "m" matrix } { "p" "positive real number" } { "n" "non-negative real number" } } +{ $description "Find the size of a matrix in 𝑙ₚ (" { $snippet "L^p" } ") vector space, usually written ∥･∥ₚ. For " { $snippet "p ≠ 1, 2, ∞" } ", this is an \"entry-wise\" norm." } +{ $examples + "Calls " { $link m-1norm } ":" + { $example + "USING: math.matrices prettyprint ;" + "4 4 1 1 matrix-p-norm ." + "4" + } + "Falls back to " { $link matrix-p-norm-entrywise } ":" + { $example + "USING: math.functions math.matrices prettyprint ;" + "2 2 3 1.5 matrix-p-norm 7.559 10e-4 ~ ." + "t" + } +} ; + +{ matrix-p-norm matrix-p-norm-entrywise } related-words +{ matrix-p-norm matrix-p-q-norm } related-words + +HELP: normalize-matrix { $values { "m" matrix } { "m'" matrix } } { $description "Normalize a matrix. Each element from the input matrix is computed as a fraction of the maximum element. The maximum element becomes " { $snippet "1/1" } "." } { $notelist @@ -948,32 +1057,11 @@ HELP: mnorm { $examples { $example "USING: math.matrices prettyprint ;" - "{ { 5 9 } { 15 17 } } mnorm ." + "{ { 5 9 } { 15 17 } } normalize-matrix ." "{ { 5/17 9/17 } { 15/17 1 } }" } } ; -HELP: m-infinity-norm -{ $values { "m" matrix } { "n" number } } ; - -HELP: m-1norm -{ $values { "m" matrix } { "n" number } } ; - -HELP: frobenius-norm -{ $values { "m" matrix } { "n" number } } -{ $notes "Also known as the Hilbert-Schmidt norm." } ; - -HELP: >square-matrix -{ $values { "m" matrix } { "subset" square-matrix } } -{ $description "Find only the " { $link2 square-matrix "square" } " subset of the input matrix." } -{ $examples - { $example - "USING: math.matrices prettyprint ;" - "{ { 0 2 4 6 } { 1 3 5 7 } } >square-matrix ." - "{ { 0 2 } { 1 3 } }" - } -} ; - HELP: main-diagonal { $values { "matrix" matrix } { "seq" sequence } } { $description "Find the main diagonal of a matrix." $nl "This diagonal begins in the upper left of the matrix at index " { $snippet "{ 0 0 }" } ", continuing downward and rightward for all indices " { $snippet "{ n n }" } " in the " { $link square-matrix } " subset of the input (see " { $link } ")." } diff --git a/basis/math/matrices/matrices-tests.factor b/basis/math/matrices/matrices-tests.factor index 7291d941c7..d55ce73616 100644 --- a/basis/math/matrices/matrices-tests.factor +++ b/basis/math/matrices/matrices-tests.factor @@ -1,4 +1,4 @@ -! Copyright (C) 2005, 2010, 2018 Slava Pestov, Joe Groff, and Cat Stevens. +! Copyright (C) 2005, 2010, 2018, 2020 Slava Pestov, Joe Groff, and Cat Stevens. USING: arrays assocs combinators.short-circuit grouping kernel math math.statistics sequences sequences.deep tools.test ; IN: math.matrices @@ -317,17 +317,6 @@ PRIVATE> mdot ] unit-test - -! TODO: note: merge conflict from HEAD contained the following -! ------------------------ -! predicates - -{ t } [ { } square-matrix? ] unit-test -{ t } [ { { 1 } } square-matrix? ] unit-test -{ t } [ { { 1 2 } { 3 4 } } square-matrix? ] unit-test -{ f } [ { { 1 } { 2 3 } } square-matrix? ] unit-test -{ f } [ { { 1 2 } } square-matrix? ] unit-test - { 9 } [ { { 2 -2 1 } { 1 3 -1 } { 2 -4 2 } } m-1norm ] unit-test @@ -336,40 +325,37 @@ PRIVATE> { 2.0 } [ { { 1 1 } { 1 1 } } frobenius-norm ] unit-test -! from "intermediate commit" -! any deep-empty matrix is null -! it doesn't make any sense for { } to be null while { { } } to be considered nonnull -{ t } [ { - { } - { { } } - { { { } } } - { { } { } { } } - { { { } } { { { } } } } -} [ null-matrix? ] all? -] unit-test -{ f } [ { - { 1 2 } - { { 1 2 } } - { { 1 } { 2 } } - { { { 1 } } { 2 } { } } -} [ null-matrix? ] any? -] unit-test +{ 10e-8 } +[ + 5.4772255 + { { 1 2 } { 3 4 } } frobenius-norm +] unit-test~ -{ t } [ 10 dup zero-matrix? ] unit-test -{ t } [ 10 10 15 zero-matrix? ] unit-test -{ t } [ 0 dup null-matrix? ] unit-test -{ f } [ 0 dup zero-matrix? ] unit-test -{ f } [ 4 zero-matrix? ] unit-test +{ 10e-6 } +[ + 36.94590 + { { 1 2 } { 4 8 } { 16 32 } } frobenius-norm +] unit-test~ -{ t } [ { } regular-matrix? ] unit-test -{ t } [ { { } } regular-matrix? ] unit-test -{ t } [ { { 1 2 } } regular-matrix? ] unit-test -{ t } [ { { 1 2 } { 3 4 } } regular-matrix? ] unit-test -{ t } [ { { 1 } { 3 } } regular-matrix? ] unit-test -{ f } [ { { 1 2 } { 3 } } regular-matrix? ] unit-test -{ f } [ { { 1 } { 3 2 } } regular-matrix? ] unit-test -! TODO: note: lines since last HEAD comment were deleted in "fix more code and add more rigorous tests" +! equivalent to frobenius for p = q = 2 +{ 2.0 } +[ { { 1 1 } { 1 1 } } 2 2 matrix-p-q-norm ] unit-test + +{ 10e-7 } +[ + 33.456466 + { { 1 2 } { 4 8 } { 16 32 } } 3 matrix-p-norm-entrywise +] unit-test~ + +{ { { -1 0 } { 0 0 } } } +[ { { -2 0 } { 0 0 } } normalize-matrix ] unit-test + +{ { { -1 0 } { 0 1/2 } } } +[ { { -2 0 } { 0 1 } } normalize-matrix ] unit-test + +{ t } +[ 3 3 dup normalize-matrix = ] unit-test ! diagonals diff --git a/basis/math/matrices/matrices.factor b/basis/math/matrices/matrices.factor index 8a64919bf1..2ca743dad0 100644 --- a/basis/math/matrices/matrices.factor +++ b/basis/math/matrices/matrices.factor @@ -1,11 +1,11 @@ -! Copyright (C) 2005, 2010, 2018 Slava Pestov, Joe Groff, and Cat Stevens. +! Copyright (C) 2005, 2010, 2018, 2020 Slava Pestov, Joe Groff, and Cat Stevens. ! See http://factorcode.org/license.txt for BSD license. USING: accessors arrays classes.singleton columns combinators combinators.short-circuit combinators.smart formatting fry grouping kernel locals math math.bits math.functions math.order math.private math.ranges math.statistics math.vectors -math.vectors.private sequences sequences.deep sequences.private -slots.private summary ; +math.vectors.private sequences sequences.deep sequences.extras +sequences.private slots.private summary ; IN: math.matrices ! defined here because of issue #1943 @@ -256,10 +256,48 @@ DEFER: matrix-set-nths : mmin ( m -- n ) [ 1/0. ] dip [ [ min ] each ] each ; : mmax ( m -- n ) [ -1/0. ] dip [ [ max ] each ] each ; -: mnorm ( m -- m' ) dup mmax abs m/n ; -: m-infinity-norm ( m -- n ) [ [ abs ] map-sum ] map supremum ; -: m-1norm ( m -- n ) flip m-infinity-norm ; -: frobenius-norm ( m -- n ) [ [ sq ] map-sum ] map-sum sqrt ; + +: m-infinity-norm ( m -- n ) + dup zero-matrix? [ drop 0. ] [ + [ [ abs ] map-sum ] map supremum + ] if ; + +: m-1norm ( m -- n ) + dup zero-matrix? [ drop 0. ] [ + flip m-infinity-norm + ] if ; + +: frobenius-norm ( m -- n ) + dup zero-matrix? [ drop 0. ] [ + [ [ sq ] map-sum ] map-sum sqrt + ] if ; + +ALIAS: hilbert-schmidt-norm frobenius-norm + +:: matrix-p-q-norm ( m p q -- n ) + m dup zero-matrix? [ drop 0. ] [ + [ [ sq ] map-sum q p / ^ ] map-sum q recip ^ + ] if ; + +: matrix-p-norm-entrywise ( m p -- n ) + dup zero-matrix? [ 2drop 0. ] [ + [ flatten1 V{ } like ] dip p-norm-default + ] if ; + +: matrix-p-norm ( m p -- n ) + dup zero-matrix? [ 2drop 0. ] [ + { + { [ dup 1 = ] [ drop m-1norm ] } + { [ dup 2 = ] [ drop frobenius-norm ] } + { [ dup fp-infinity? ] [ drop m-infinity-norm ] } + [ matrix-p-norm-entrywise ] + } cond + ] if ; + +: normalize-matrix ( m -- m' ) + dup zero-matrix? [ ] [ + dup mabs mmax m/n + ] if ; ! well-defined for square matrices; but works on nonsquare too : main-diagonal ( matrix -- seq ) diff --git a/extra/math/matrices/extras/extras-docs.factor b/extra/math/matrices/extras/extras-docs.factor index b16f015a0a..2a14d895f1 100644 --- a/extra/math/matrices/extras/extras-docs.factor +++ b/extra/math/matrices/extras/extras-docs.factor @@ -49,9 +49,6 @@ ARTICLE: "math.matrices.extras" "Extra matrix operations" "Matrix algebra:" { $subsections - mmin - mmax - mnorm rank nullity