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.<hilbert-matrix>`
	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.
master
Cat Stevens 2020-02-29 21:46:38 -05:00 committed by John Benediktsson
parent 2511fa72de
commit 573e4ed198
4 changed files with 186 additions and 77 deletions

View File

@ -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 <matrix> 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 <matrix> 1 matrix-p-norm ."
"4"
}
"Falls back to " { $link matrix-p-norm-entrywise } ":"
{ $example
"USING: math.functions math.matrices prettyprint ;"
"2 2 3 <matrix> 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 <square-rows> } ")." }

View File

@ -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> zero-matrix? ] unit-test
{ t } [ 10 10 15 <simple-eye> zero-matrix? ] unit-test
{ t } [ 0 dup <zero-matrix> null-matrix? ] unit-test
{ f } [ 0 dup <zero-matrix> zero-matrix? ] unit-test
{ f } [ 4 <identity-matrix> 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 <zero-matrix> dup normalize-matrix = ] unit-test
! diagonals

View File

@ -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 )

View File

@ -49,9 +49,6 @@ ARTICLE: "math.matrices.extras" "Extra matrix operations"
"Matrix algebra:"
{ $subsections
mmin
mmax
mnorm
rank
nullity