Opened 7 years ago
Closed 3 years ago
#13703 closed enhancement (fixed)
special matrices
Reported by:  jason  Owned by:  jason, was 

Priority:  minor  Milestone:  sage8.0 
Component:  linear algebra  Keywords:  
Cc:  rbeezer, kcrisman, mforets  Merged in:  
Authors:  Jason Grout, Marcelo Forets  Reviewers:  Dima Pasechnik, Marcelo Forets 
Report Upstream:  N/A  Work issues:  
Branch:  b29d152 (Commits)  Commit:  b29d152bdafed69f6030a1870f38587177995d84 
Dependencies:  Stopgaps: 
Description (last modified by )
It would be great to have a matrices namespace to put special matrix commands, sort of like the:
 graphs.* namespace
 groups.* namespace
 polytopes.* namespace
Here are some starter definitions:
def hilbert(R,n): return matrix(R, n, lambda i,j: 1/(i+j+1)) def vandermonde(R, v): return matrix(R, len(v), lambda i,j: v[i]^j) def toeplitz(R,c,r): return matrix(R, len(c), len(r), lambda i,j: c[ij] if i>=j else r[ji]) def hankel(R,c,r): entries=c+r[1:]; return matrix(R, len(c), len(r), lambda i,j: entries[i+j]) def circulant(R, E): return toeplitz(R, E[0:1]+E[1:0:1], E) def skew_circulant(R,E): return hankel(R, E, E[1:]+E[:1]) #Hadamard matrices: def legendre_symbol(x): """Extend the built in legendre_symbol function to handle prime power fields. Assume x is an element of a finite field as well""" if x==0: return 0 elif x.is_square(): return 1 else: return 1 def jacobsthal(p,n): """See http://en.wikipedia.org/wiki/Paley_construction for a way to use jacobsthal matrices to construct hadamard matrices""" if n == 1: elts = GF(p).list() else: elts = GF(p^n,'a').list() return matrix(len(elts), lambda i,j: legendre_symbol(elts[i]elts[j])) def paley_matrix(p,n): """See http://en.wikipedia.org/wiki/Paley_construction""" mod = p^n%4 if mod == 3: # Paley Type 1 construction ones = vector([1]*p^n) QplusI = jacobsthal(p,n) # Q+=I efficiently for i in range(p^n): QplusI[i,i]=1 return block_matrix(2,[ [1,ones.row()], [ones.column(), QplusI]]) elif mod == 1: # Paley Type 2 construction ones = vector([1]*p^n) QplusI = jacobsthal(p,n) QminusI = copy(QplusI) for i in range(p^n): QplusI[i,i]=1 QminusI[i,i]=1 SplusI = block_matrix(2,[[1,ones.row()],[ones.column(), QplusI]]) SminusI = block_matrix(2,[[1,ones.row()], [ones.column(), QminusI]]) return block_matrix(2,[[SplusI,SminusI],[SminusI,SplusI]]) else: raise ValueError("p^n must be congruent to 1 or 3 mod 4") }} The matrixform patch below first implements a decorator to take care of an optional first argument of the ring, as well as an optional keyword argument of 'sparse'. Then it implements most of these using that decorator, so the boilerplate code is abstracted away. Additionally, we could use scipy to create more matrices (or do it ourselves): http://docs.scipy.org/doc/scipy/reference/linalg.html#specialmatrices (thanks to pascal on sagesupport for correcting the circulant code above: https://groups.google.com/d/msg/sagesupport/RnKjQ9n2YB0/vfCEvIV_HZUJ )
Attachments (2)
Change History (57)
comment:1 Changed 7 years ago by
comment:2 Changed 7 years ago by
Ah, yes, nice. Should we close this ticket and move these examples over or make it an addition to that ticket?
comment:3 Changed 7 years ago by
I'd vote for implementing the new matrices on this ticket, making sure they build on the namespace/tabcompletion decorator in #13678. Just change the description slightly to speak of just the special matrices.
Note: These could be named, e.g., hilbert_matrix()
, and the _matrix
part gets stripped put when found with tabcompletion.
comment:4 Changed 7 years ago by
 Description modified (diff)
comment:5 Changed 7 years ago by
 Description modified (diff)
More code, this time for Hadamard matrices. Thanks to wikipedia and a paper by Alex Kramer summarizing these ideas. I can't seem to find a source for the Paley type 2 matrix construction, but it seems to check out for small values of p^{n, and it looks like a reasonable modification of the wikipedia method. }
comment:6 Changed 7 years ago by
Maybe just convert the example definitions to actual functions, add a few examples, and add them to Sage. It's better than letting this code "bitrot" in trac :)
comment:7 Changed 7 years ago by
 Description modified (diff)
Changed 7 years ago by
comment:8 Changed 7 years ago by
 Status changed from new to needs_review
comment:9 Changed 7 years ago by
 Status changed from needs_review to needs_work
 Work issues set to documentation, doctests, checking
Documentation (and doctests) should still be added to the new functions, as well as the offset argument for the diagonal matrix constructor.
Ideally, we should have a new decorator that handles the optional ring argument automatically.
comment:10 Changed 7 years ago by
 Cc kcrisman added
Adding that documentation sounds like a great student project...
comment:11 Changed 7 years ago by
 Description modified (diff)
comment:12 Changed 7 years ago by
I would be happy to review this patch. Should I ?
comment:13 Changed 7 years ago by
The functions still need documentation, so it's not quite ready to go in yet. However, if you want to add docs, or if you want to check correctness of the functions, that would be great too.
comment:14 Changed 6 years ago by
 Milestone changed from sage5.11 to sage5.12
comment:15 Changed 6 years ago by
 Milestone changed from sage6.1 to sage6.2
comment:16 Changed 6 years ago by
 Milestone changed from sage6.2 to sage6.3
comment:17 Changed 5 years ago by
 Milestone changed from sage6.3 to sage6.4
comment:18 Changed 3 years ago by
 Cc mforets added
 Description modified (diff)
comment:19 Changed 3 years ago by
 Branch set to u/mforets/13703
 Commit set to 039b3e29cff022f820a5299479f5f4f4db5f9ed4
 Milestone changed from sage6.4 to sage8.0
 added documentation and tests to hankel, vandermonde, toeplitz and hilbert
 using the approach from the 1st patch (matrix method) decorator, because it is more in line with the current
special.py
 added the file
special.py
to the html reference, because it was not there at all  added a table with the available constructions, in the same lines as with the polytopes library
it is missing to document and add doctests for the Paley construction. perhaps a person with a number theory background can continue that part.

this is related to documentation here but less important issue: since we are using the matrix method decorator, it adds the strings This function is available as xxx(...) and matrix.xxx(...).
, although this is not true in general. it applies not only for this ticket; see for instance the help for random_rref_matrix
.
New commits:
8ccddab  add doc and examples to 4 functions

039b3e2  switch to matix method approach

comment:20 followup: ↓ 23 Changed 3 years ago by
note that there is now a lot on Hadamard matrices in src/sage/combinat/matrices/hadamard_matrix.py
so no need to duplicate this in special.py
. (Perhaps more docs can be added and/or crosslinked)
comment:21 followup: ↓ 22 Changed 3 years ago by
regarding is available as xxx(...) and matrix.xxx(...)
, I think there is a consensus that nothing should be added to the global namespace, which is already quite full, rather, it should be pruned...
comment:22 in reply to: ↑ 21 Changed 3 years ago by
Replying to dimpase:
regarding
is available as xxx(...) and matrix.xxx(...)
, I think there is a consensus that nothing should be added to the global namespace, which is already quite full, rather, it should be pruned...
the @matrix_method
decorator adds the possibility for tabcompletion, hence matrix.xxx
works. this is good. but it also adds a string to the inline help saying that
... prefix = " This function is available as %s(...) and matrix.%s(...)." % (func.__name__, name ...
but the %s(...)
is not true; it is only true if the function is included in all.py
.
i was not proposing to add anything to all.py
, just that someone revises and eventually confirms that this prefix to the docstring is false and should be fixed (or not).
comment:23 in reply to: ↑ 20 Changed 3 years ago by
Replying to dimpase:
note that there is now a lot on Hadamard matrices in
src/sage/combinat/matrices/hadamard_matrix.py
so no need to duplicate this inspecial.py
. (Perhaps more docs can be added and/or crosslinked)
this is great, thanks for pointing it out! i'll take care to add the relevant crosslink, in the new class docstring of special.py
comment:24 Changed 3 years ago by
 Commit changed from 039b3e29cff022f820a5299479f5f4f4db5f9ed4 to 84d105685e556a05540f11ba91589e40c0a242ac
comment:25 Changed 3 years ago by
 Status changed from needs_work to needs_review
 Work issues documentation, doctests, checking deleted
 added Jason Grout as author, because this applies his patch to 4 special matrices (uploaded the original patches)
 leaving the (minor) issue that i didn't understand about the meaning of the
prefix
string of the@matrix_method
decorator outside the scope of this ticket.  adding myself as reviewer (?)
 added cross links to combinat
can someone else have a look? thanks !
comment:26 Changed 3 years ago by
 Commit changed from 84d105685e556a05540f11ba91589e40c0a242ac to 4fdccb6cf9a33d136521b5d6ada06b43eae9f7f1
Branch pushed to git repo; I updated commit sha1. New commits:
4fdccb6  revert to original module title

comment:27 Changed 3 years ago by
 Reviewers set to Dima Pasechnik, Marcelo Forets
comment:28 Changed 3 years ago by
 Status changed from needs_review to positive_review
looks good to me.
comment:29 Changed 3 years ago by
 Status changed from positive_review to needs_work
Thanks dimpase for the review!!! I've realized the word "vector" is innacurate in the docstring of hankel, because entries = c + r[1:]
makes sense for lists.
do you have some recommendation about being permissible with arguments? (only list/only vector/both). from my side i'll check how other parts of sage handle a similar issue.
comment:30 Changed 3 years ago by
 Commit changed from 4fdccb6cf9a33d136521b5d6ada06b43eae9f7f1 to 4d8e6bb2daa1784a8f5c203bd596e80b6e05e8aa
Branch pushed to git repo; I updated commit sha1. New commits:
4d8e6bb  accept vector in hankel, typos

comment:31 Changed 3 years ago by
 Status changed from needs_work to needs_review
propose to do the computation by accessing directly to the elements (as in the other cases). this way without type checking we are a little more robust in the sense that it accepts both lists and vectors. it is mathematically equivalent.
there is a price to pay indeed if we use a lambda function instead of straightforward concatenation of lists, but i judge "small" (~30ms for a 300x600 hankel matrix).
in sage there are indeed a couple of docstrings which explicitly say "vector or list" (or "list or vector") as a search shows, eg. setting rows/cols in a matrix, or a polyhedral fan. in our case i'm just keeping the word "vector"..
comment:32 followup: ↓ 33 Changed 3 years ago by
actually, I don't like the interface for Hankel matrices. The normal convention is that you have a potentially infinite matrix
a b c d e f..... b c d e f ...... c d e f ........ ................
corresponding to the sequence (a,b,c,d,e,f,...) and then you just cut out the first n rows and columns. (so it's square by default), see https://en.wikipedia.org/wiki/Hankel_matrix
comment:33 in reply to: ↑ 32 Changed 3 years ago by
Replying to dimpase:
actually, I don't like the interface for Hankel matrices. The normal convention is that you have a potentially infinite matrix
a b c d e f..... b c d e f ...... c d e f ........ ................corresponding to the sequence (a,b,c,d,e,f,...) and then you just cut out the first n rows and columns. (so it's square by default), see https://en.wikipedia.org/wiki/Hankel_matrix
ok. square by default is:
def hankel(c, r=None): m = len(c) if r is not None: n = len(r) entries = lambda i: c[i] if i < m else r[im+1] else: n = len(c) entries = lambda i: c[i] if i < m else c[im+1] return matrix(lambda i,j: entries(i+j), nrows=m, ncols=n)
i don't have any particular opinion about allowing rectangular matrices as optional, or removing it altogether.
notice that the original code was like that, and there a similar Matlab command: hankel(c, r)
.
comment:34 Changed 3 years ago by
 Commit changed from 4d8e6bb2daa1784a8f5c203bd596e80b6e05e8aa to c636bbfe1576b124b15e795b51570e5a1e382fca
Branch pushed to git repo; I updated commit sha1. New commits:
c636bbf  square by default in Hankel

comment:35 Changed 3 years ago by
 Commit changed from c636bbfe1576b124b15e795b51570e5a1e382fca to a5c5ed68440581ae7f4df1b8f5f3abe7b82073e6
Branch pushed to git repo; I updated commit sha1. New commits:
a5c5ed6  fix typo

comment:36 Changed 3 years ago by
@dimpase : done with the suggested interface to the Hankel matrices. let me know if you'd like me to add other special matrices or related issues to this ticket (eg. custom specification of the base_ring
).
something i have needed before (and there are several requests for this in asksage as you know) is to provide an "abstract matrix" constructor like matrix.symbolic(5)
or matrix.abstract('a', 5)
to return the matrix with SR coefficients labelled aij
.
comment:37 Changed 3 years ago by
 Commit changed from a5c5ed68440581ae7f4df1b8f5f3abe7b82073e6 to b58e47e1f1fefc0093d831a1580beda2bc9ad424
comment:38 Changed 3 years ago by
from my side this ticket is good to go :) i've run tests and built the html doc. please let me now if you have changes or improvements so i can make them. Thanks.
comment:39 followup: ↓ 41 Changed 3 years ago by
Please add an example illustrating the rule about row/column "conflict" for Hankel matrices.
comment:40 Changed 3 years ago by
 Commit changed from b58e47e1f1fefc0093d831a1580beda2bc9ad424 to 5f4504bb89d8f75321d0caa5c6fa38fa41ab5913
Branch pushed to git repo; I updated commit sha1. New commits:
5f4504b  add two illustrating examples in Hankel

comment:41 in reply to: ↑ 39 Changed 3 years ago by
Replying to dimpase:
Please add an example illustrating the rule about row/column "conflict" for Hankel matrices.
Good idea, done.
comment:42 followup: ↓ 43 Changed 3 years ago by
OK, at least now I understand the rule implemented easily. But why is this chosen?
Would it be more natural to just ignore the "conflict" thing completely and assume that the the 2nd argument is the last row without the 1st entry? This way the users do not even need to know the last entry of the 1st column to have full control over the matrix they create. E.g., in the present design one can have a match by accident, resulting in an wrong matrix. In the design I propose the user does not need to worry about this.
comment:43 in reply to: ↑ 42 Changed 3 years ago by
Replying to dimpase:
OK, at least now I understand the rule implemented easily. But why is this chosen?
it was there in the patch in this thread, and it seemed ok since it is similar in matlab software, scipy, and mathematica. But ok to improve things.
Would it be more natural to just ignore the "conflict" thing completely and assume that the the 2nd argument is the last row without the 1st entry? This way the users do not even need to know the last entry of the 1st column to have full control over the matrix they create. E.g., in the present design one can have a match by accident, resulting in an wrong matrix. In the design I propose the user does not need to worry about this.
sorry i didnt't follow the suggestion; do you mean at the docstring level?
(e.g. conflict > ignore
? that sounds good) if it's in the code could you reformulate please? maybe an example would help.
comment:44 followup: ↓ 47 Changed 3 years ago by
In the branch implementation ofhankel(a,b)
, the 1st entry of b
, i.e. b[0]
is always ignored. While this is what Matlab does (Mathematica too, and Mathematica documentation requires the last entry of a
to be equal to the 1st entry of b
), it is ugly and illogical.
I propose that hankel(a,b)
should do the same what hankel(a,[None]+b)
does in the implementation on the branch.
Same applies to toeplitz()
 it suffers from the same inconsistency.
comment:45 Changed 3 years ago by
ok, we are on the same page now. i'll take care.
comment:46 Changed 3 years ago by
 Commit changed from 5f4504bb89d8f75321d0caa5c6fa38fa41ab5913 to 9d92967f4eca938b8446ca21a0fbd0f4bfeb9621
Branch pushed to git repo; I updated commit sha1. New commits:
f39ec3c  change hankel interface

3824638  add base ring to hankel

5950b5e  add base ring to INPUT of other functions

0b1da36  Merge branch 'develop' into t/mforets/13703/matrixforms

c3a8a67  change Toeplitz design; add one more example

9d92967  precision in toeplitz docstring

comment:47 in reply to: ↑ 44 Changed 3 years ago by
Replying to dimpase:
In the branch implementation of
hankel(a,b)
, the 1st entry ofb
, i.e.b[0]
is always ignored. While this is what Matlab does (Mathematica too, and Mathematica documentation requires the last entry ofa
to be equal to the 1st entry ofb
), it is ugly and illogical.I propose that
hankel(a,b)
should do the same whathankel(a,[None]+b)
does in the implementation on the branch.Same applies to
toeplitz()
 it suffers from the same inconsistency.
done ==> needs review
comment:48 followup: ↓ 50 Changed 3 years ago by
this is nitpicking, I know, sorry, but why do you prepend [None] to r
in
+ r = c if r is None else [None] + list(r) + n = len(r) + entries = lambda i: c[i] if i < m else r[im+1]
instead of adjusting the indexing (i.e. removing +1
) of r
in the last line (and adding 1
to n
)?
comment:49 Changed 3 years ago by
While we are at it, I don't understand the choice made for r if it is not supplied. Would all 0s be more natural, i.e.? (This is what MMa's do, in fact).
if r is None: r = [0*c[0]]*(m1)
comment:50 in reply to: ↑ 48 Changed 3 years ago by
no worries, i'm here to learn.
if instead of the code in the branch one has:
sage: r = c if r is None else list(r) sage: n = len(r) sage: entries = lambda i: c[i] if i < m else r[im] sage: return matrix(lambda i,j: entries(i+j), nrows=m, ncols=n+1, ring=ring)
then the behaviour when r is not given does change (of course the matrix is no longer square):
Failed example: matrix.hankel([1..3]) Expected: [1 2 3] [2 3 2] [3 2 3] Got: [1 2 3 1] [2 3 1 2] [3 1 2 3]
i think this was the reason i kept the None.
however and in relation to your newer comment, i understand correctly, you've suggested to have this:
sage: r = c[1:] if r is None else list(r) sage: n = len(r) sage: entries = lambda i: c[i] if i < m else r[im] sage: return matrix(lambda i,j: entries(i+j), nrows=m, ncols=n+1, ring=ring)
in this situation, all current tests pass.
comment:51 Changed 3 years ago by
well, sorry i read too fast :/ let me think about all zeros.
comment:52 Changed 3 years ago by
 Commit changed from 9d92967f4eca938b8446ca21a0fbd0f4bfeb9621 to b29d152bdafed69f6030a1870f38587177995d84
Branch pushed to git repo; I updated commit sha1. New commits:
b29d152  change hankel 0's below and update doctests

comment:53 Changed 3 years ago by
To summarize:
the current behaviour for hankel(c)
is the same as in matlab and mathematica. for example,
sage: matrix.hankel([1, 2, 3, 4]) [1 2 3 4] [2 3 4 0] [3 4 0 0] [4 0 0 0] sage: matlab('hankel([1, 2, 3, 4])') 1 2 3 4 2 3 4 0 3 4 0 0 4 0 0 0
however, if row and column are passed all the data in the row is considered:
sage: matrix.hankel([1, 2, 3], [7, 8, 9, 10]) [ 1 2 3 7 8] [ 2 3 7 8 9] [ 3 7 8 9 10]
while:
sage: matlab('hankel([1, 2, 3], [7, 8, 9, 10])') 1 2 3 8 2 3 8 9 3 8 9 10
comment:54 Changed 3 years ago by
 Status changed from needs_review to positive_review
Looks good to me.
comment:55 Changed 3 years ago by
 Branch changed from u/mforets/13703 to b29d152bdafed69f6030a1870f38587177995d84
 Resolution set to fixed
 Status changed from positive_review to closed
Like this: #13678 ? (Which I should finish reviewing.)