Opened 7 years ago

Closed 3 years ago

#13703 closed enhancement (fixed)

special matrices

Reported by: jason Owned by: jason, was
Priority: minor Milestone: sage-8.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 mforets)

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[i-j] if i>=j else r[j-i])
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#special-matrices

(thanks to pascal on sage-support for correcting the circulant code above: https://groups.google.com/d/msg/sage-support/RnKjQ9n2YB0/vfCEvIV_HZUJ )

Attachments (2)

trac-13703.patch (4.4 KB) - added by jason 7 years ago.
trac-13703-matrixform.patch (5.2 KB) - added by jason 7 years ago.
take 2, apply this patch *or* the other patch, but not both.

Download all attachments as: .zip

Change History (57)

comment:1 Changed 7 years ago by rbeezer

Like this: #13678 ? (Which I should finish reviewing.)

comment:2 Changed 7 years ago by jason

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 rbeezer

I'd vote for implementing the new matrices on this ticket, making sure they build on the namespace/tab-completion 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 tab-completion.

comment:4 Changed 7 years ago by jason

  • Description modified (diff)

comment:5 Changed 7 years ago by jason

  • 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 pn, and it looks like a reasonable modification of the wikipedia method.

comment:6 Changed 7 years ago by ppurka

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 jason

  • Description modified (diff)

Changed 7 years ago by jason

comment:8 Changed 7 years ago by jason

  • Status changed from new to needs_review

comment:9 Changed 7 years ago by jason

  • 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 kcrisman

  • Cc kcrisman added

Adding that documentation sounds like a great student project...

Changed 7 years ago by jason

take 2, apply this patch *or* the other patch, but not both.

comment:11 Changed 7 years ago by jason

  • Description modified (diff)

comment:12 Changed 7 years ago by Bouillaguet

I would be happy to review this patch. Should I ?

comment:13 Changed 7 years ago by jason

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 jdemeyer

  • Milestone changed from sage-5.11 to sage-5.12

comment:15 Changed 6 years ago by vbraun_spam

  • Milestone changed from sage-6.1 to sage-6.2

comment:16 Changed 6 years ago by vbraun_spam

  • Milestone changed from sage-6.2 to sage-6.3

comment:17 Changed 5 years ago by vbraun_spam

  • Milestone changed from sage-6.3 to sage-6.4

comment:18 Changed 3 years ago by mforets

  • Cc mforets added
  • Description modified (diff)

comment:19 Changed 3 years ago by mforets

  • Branch set to u/mforets/13703
  • Commit set to 039b3e29cff022f820a5299479f5f4f4db5f9ed4
  • Milestone changed from sage-6.4 to sage-8.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:

8ccddabadd doc and examples to 4 functions
039b3e2switch to matix method approach

comment:20 follow-up: Changed 3 years ago by 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 in special.py. (Perhaps more docs can be added and/or cross-linked)

comment:21 follow-up: Changed 3 years ago by 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...

comment:22 in reply to: ↑ 21 Changed 3 years ago by mforets

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 tab-completion, 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 mforets

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 in special.py. (Perhaps more docs can be added and/or cross-linked)

this is great, thanks for pointing it out! i'll take care to add the relevant cross-link, in the new class docstring of special.py

comment:24 Changed 3 years ago by git

  • Commit changed from 039b3e29cff022f820a5299479f5f4f4db5f9ed4 to 84d105685e556a05540f11ba91589e40c0a242ac

Branch pushed to git repo; I updated commit sha1. New commits:

a08f04aremove unused decorator
84d1056add cross-link to combinats module

comment:25 Changed 3 years ago by mforets

  • Authors set to Jason Grout
  • 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 !

Last edited 3 years ago by mforets (previous) (diff)

comment:26 Changed 3 years ago by git

  • Commit changed from 84d105685e556a05540f11ba91589e40c0a242ac to 4fdccb6cf9a33d136521b5d6ada06b43eae9f7f1

Branch pushed to git repo; I updated commit sha1. New commits:

4fdccb6revert to original module title

comment:27 Changed 3 years ago by dimpase

  • Authors changed from Jason Grout to Jason Grout, Marcelo Forets
  • Reviewers set to Dima Pasechnik, Marcelo Forets

comment:28 Changed 3 years ago by dimpase

  • Status changed from needs_review to positive_review

looks good to me.

comment:29 Changed 3 years ago by mforets

  • 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 git

  • Commit changed from 4fdccb6cf9a33d136521b5d6ada06b43eae9f7f1 to 4d8e6bb2daa1784a8f5c203bd596e80b6e05e8aa

Branch pushed to git repo; I updated commit sha1. New commits:

4d8e6bbaccept vector in hankel, typos

comment:31 Changed 3 years ago by mforets

  • 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 follow-up: Changed 3 years ago by 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

comment:33 in reply to: ↑ 32 Changed 3 years ago by mforets

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[i-m+1]
    else:
        n = len(c)
        entries = lambda i: c[i] if i < m else c[i-m+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 git

  • Commit changed from 4d8e6bb2daa1784a8f5c203bd596e80b6e05e8aa to c636bbfe1576b124b15e795b51570e5a1e382fca

Branch pushed to git repo; I updated commit sha1. New commits:

c636bbfsquare by default in Hankel

comment:35 Changed 3 years ago by git

  • Commit changed from c636bbfe1576b124b15e795b51570e5a1e382fca to a5c5ed68440581ae7f4df1b8f5f3abe7b82073e6

Branch pushed to git repo; I updated commit sha1. New commits:

a5c5ed6fix typo

comment:36 Changed 3 years ago by mforets

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

  • Commit changed from a5c5ed68440581ae7f4df1b8f5f3abe7b82073e6 to b58e47e1f1fefc0093d831a1580beda2bc9ad424

Branch pushed to git repo; I updated commit sha1. New commits:

c978b57use optional ring keyword
c847840Merge branch 'develop' into t/mforets/13703/matrixforms
b58e47esome tweaks

comment:38 Changed 3 years ago by mforets

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 follow-up: Changed 3 years ago by dimpase

Please add an example illustrating the rule about row/column "conflict" for Hankel matrices.

comment:40 Changed 3 years ago by git

  • Commit changed from b58e47e1f1fefc0093d831a1580beda2bc9ad424 to 5f4504bb89d8f75321d0caa5c6fa38fa41ab5913

Branch pushed to git repo; I updated commit sha1. New commits:

5f4504badd two illustrating examples in Hankel

comment:41 in reply to: ↑ 39 Changed 3 years ago by mforets

Replying to dimpase:

Please add an example illustrating the rule about row/column "conflict" for Hankel matrices.

Good idea, done.

comment:42 follow-up: Changed 3 years ago by dimpase

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 mforets

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 follow-up: Changed 3 years ago by dimpase

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 mforets

ok, we are on the same page now. i'll take care.

comment:46 Changed 3 years ago by git

  • Commit changed from 5f4504bb89d8f75321d0caa5c6fa38fa41ab5913 to 9d92967f4eca938b8446ca21a0fbd0f4bfeb9621

Branch pushed to git repo; I updated commit sha1. New commits:

f39ec3cchange hankel interface
3824638add base ring to hankel
5950b5eadd base ring to INPUT of other functions
0b1da36Merge branch 'develop' into t/mforets/13703/matrixforms
c3a8a67change Toeplitz design; add one more example
9d92967precision in toeplitz docstring

comment:47 in reply to: ↑ 44 Changed 3 years ago by mforets

Replying to dimpase:

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.

done ==> needs review

comment:48 follow-up: Changed 3 years ago by dimpase

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[i-m+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 dimpase

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]]*(m-1)
Last edited 3 years ago by dimpase (previous) (diff)

comment:50 in reply to: ↑ 48 Changed 3 years ago by mforets

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[i-m]
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[i-m]
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 mforets

well, sorry i read too fast :/ let me think about all zeros.

comment:52 Changed 3 years ago by git

  • Commit changed from 9d92967f4eca938b8446ca21a0fbd0f4bfeb9621 to b29d152bdafed69f6030a1870f38587177995d84

Branch pushed to git repo; I updated commit sha1. New commits:

b29d152change hankel 0's below and update doctests

comment:53 Changed 3 years ago by mforets

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 dimpase

  • Status changed from needs_review to positive_review

Looks good to me.

comment:55 Changed 3 years ago by vbraun

  • Branch changed from u/mforets/13703 to b29d152bdafed69f6030a1870f38587177995d84
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.