Opened 9 years ago

Closed 5 years ago

# special matrices

Reported by: Owned by: jason jason, was minor sage-8.0 linear algebra rbeezer, kcrisman, mforets Jason Grout, Marcelo Forets Dima Pasechnik, Marcelo Forets N/A b29d152 b29d152bdafed69f6030a1870f38587177995d84

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 skew_circulant(R,E): return hankel(R, E, E[-1:]+E[:-1])

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

### comment:1 Changed 9 years ago by rbeezer

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

### comment:2 Changed 9 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 9 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 9 years ago by jason

• Description modified (diff)

### comment:5 Changed 9 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 9 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 9 years ago by jason

• Description modified (diff)

### comment:8 Changed 9 years ago by jason

• Status changed from new to needs_review

### comment:9 Changed 9 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 9 years ago by kcrisman

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

### Changed 9 years ago by jason

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

### comment:11 Changed 9 years ago by jason

• Description modified (diff)

### comment:12 Changed 9 years ago by Bouillaguet

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

### comment:13 Changed 9 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 8 years ago by jdemeyer

• Milestone changed from sage-5.11 to sage-5.12

### comment:15 Changed 8 years ago by vbraun_spam

• Milestone changed from sage-6.1 to sage-6.2

### comment:16 Changed 8 years ago by vbraun_spam

• Milestone changed from sage-6.2 to sage-6.3

### comment:17 Changed 7 years ago by vbraun_spam

• Milestone changed from sage-6.3 to sage-6.4

### comment:18 Changed 5 years ago by mforets

• Description modified (diff)

### comment:19 Changed 5 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:

 ​8ccddab `add doc and examples to 4 functions` ​039b3e2 `switch to matix method approach`

### comment:20 follow-up: ↓ 23 Changed 5 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: ↓ 22 Changed 5 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 5 years ago by mforets

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 5 years ago by mforets

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 5 years ago by git

• Commit changed from 039b3e29cff022f820a5299479f5f4f4db5f9ed4 to 84d105685e556a05540f11ba91589e40c0a242ac

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

 ​a08f04a `remove unused decorator` ​84d1056 `add cross-link to combinats module`

### comment:25 Changed 5 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 (?)

can someone else have a look? thanks !

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

### comment:26 Changed 5 years ago by git

• 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 5 years ago by dimpase

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

### comment:28 Changed 5 years ago by dimpase

• Status changed from needs_review to positive_review

looks good to me.

### comment:29 Changed 5 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 5 years ago by git

• 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 5 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: ↓ 33 Changed 5 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 5 years ago by mforets

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 5 years ago by git

• 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 5 years ago by git

• Commit changed from c636bbfe1576b124b15e795b51570e5a1e382fca to a5c5ed68440581ae7f4df1b8f5f3abe7b82073e6

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

 ​a5c5ed6 `fix typo`

### comment:36 Changed 5 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 5 years ago by git

• Commit changed from a5c5ed68440581ae7f4df1b8f5f3abe7b82073e6 to b58e47e1f1fefc0093d831a1580beda2bc9ad424

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

 ​c978b57 `use optional ring keyword` ​c847840 `Merge branch 'develop' into t/mforets/13703/matrixforms` ​b58e47e `some tweaks`

### comment:38 Changed 5 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:40 Changed 5 years ago by git

• Commit changed from b58e47e1f1fefc0093d831a1580beda2bc9ad424 to 5f4504bb89d8f75321d0caa5c6fa38fa41ab5913

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

 ​5f4504b `add two illustrating examples in Hankel`

Good idea, done.

### comment:42 follow-up: ↓ 43 Changed 5 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 5 years ago by mforets

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: ↓ 47 Changed 5 years ago by dimpase

In the branch implementation of`hankel(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 5 years ago by mforets

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

### comment:46 Changed 5 years ago by git

• 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 5 years ago by mforets

In the branch implementation of`hankel(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: ↓ 50 Changed 5 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 5 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 5 years ago by dimpase (previous) (diff)

### comment:50 in reply to: ↑ 48 Changed 5 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 5 years ago by mforets

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

### comment:52 Changed 5 years ago by git

• 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 5 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 5 years ago by dimpase

• Status changed from needs_review to positive_review

Looks good to me.

### comment:55 Changed 5 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.