linalg.toeplitz and ticker #667
I am going to update linalg.toeplitz to fix the problem reported in ticket #667, but I have some questions about the desired behavior of the function when either `c` (the first column) or `r` (the first row) is complex. The docstring does not say anything about complex values, but the code does some undocumented and unexpected things in this case, one example of which is the bug reported in the ticket #667. I would like to fix this, but I would first like opinions on exactly how it *should* handle complex arguments. I *think* the idea is that if `c` is complex and `r` is not given, then a Hermitian Toeplitz matrix is created. Is that correct? Note that this requires that c[0] be real, since c[0] is the value that will be in the diagonal of the matrix. While I am fixing this, I would like to remove the printed warning that occurs when r[0] != c[0]. I would like to change this to raise a ValueError--that is, adopt the philosophy that either the arguments are correct (and the code handles them as [not yet, but soon to be] documented in the docstring, with no warnings needed), or the arguments are wrong and an exception should be raised. This means that when `r` is not given, an exception would be raised if c[0] is not real. A very different way to handle the case of `r` not being given is to assume it means a square matrix should be created with r[1:] = 0 (i.e. zeros in all the upper off-diagonals). This avoids the whole issue of complex numbers and conjugates, but it is also a more drastic change. Warren
Warren Weckesser wrote:
I am going to update linalg.toeplitz to fix the problem reported in ticket #667, but I have some questions about the desired behavior of the function when either `c` (the first column) or `r` (the first row) is complex. The docstring does not say anything about complex values, but the code does some undocumented and unexpected things in this case, one example of which is the bug reported in the ticket #667.
I would like to fix this, but I would first like opinions on exactly how it *should* handle complex arguments. I *think* the idea is that if `c` is complex and `r` is not given, then a Hermitian Toeplitz matrix is created. Is that correct? Note that this requires that c[0] be real, since c[0] is the value that will be in the diagonal of the matrix.
While I am fixing this, I would like to remove the printed warning that occurs when r[0] != c[0]. I would like to change this to raise a ValueError--that is, adopt the philosophy that either the arguments are correct (and the code handles them as [not yet, but soon to be] documented in the docstring, with no warnings needed), or the arguments are wrong and an exception should be raised.
This means that when `r` is not given, an exception would be raised if c[0] is not real.
A very different way to handle the case of `r` not being given is to assume it means a square matrix should be created with r[1:] = 0 (i.e. zeros in all the upper off-diagonals). This avoids the whole issue of complex numbers and conjugates, but it is also a more drastic change.
By the way, that is how linalg.hankel handles the case where `r` is not given. Warren
On Sun, Mar 7, 2010 at 8:32 PM, Warren Weckesser <warren.weckesser@enthought.com> wrote:
Warren Weckesser wrote:
I am going to update linalg.toeplitz to fix the problem reported in ticket #667, but I have some questions about the desired behavior of the function when either `c` (the first column) or `r` (the first row) is complex. The docstring does not say anything about complex values, but the code does some undocumented and unexpected things in this case, one example of which is the bug reported in the ticket #667.
I would like to fix this, but I would first like opinions on exactly how it *should* handle complex arguments. I *think* the idea is that if `c` is complex and `r` is not given, then a Hermitian Toeplitz matrix is created. Is that correct? Note that this requires that c[0] be real, since c[0] is the value that will be in the diagonal of the matrix.
While I am fixing this, I would like to remove the printed warning that occurs when r[0] != c[0]. I would like to change this to raise a ValueError--that is, adopt the philosophy that either the arguments are correct (and the code handles them as [not yet, but soon to be] documented in the docstring, with no warnings needed), or the arguments are wrong and an exception should be raised.
This means that when `r` is not given, an exception would be raised if c[0] is not real.
A very different way to handle the case of `r` not being given is to assume it means a square matrix should be created with r[1:] = 0 (i.e. zeros in all the upper off-diagonals). This avoids the whole issue of complex numbers and conjugates, but it is also a more drastic change.
a lot of code relies on toeplitz of a single argument to build the symmetric matrix.
By the way, that is how linalg.hankel handles the case where `r` is not given.
??
import scipy.linalg scipy.linalg.toeplitz(np.arange(4)) array([[0, 1, 2, 3], [1, 0, 1, 2], [2, 1, 0, 1], [3, 2, 1, 0]]) scipy.linalg.hankel(np.arange(4)) array([[ 0., 1., 2., 3.], [ 1., 2., 3., 0.], [ 2., 3., 0., 0.], [ 3., 0., 0., 0.]])
scipy.linalg.hankel(np.arange(4)+1.j) array([[ 0.+1.j, 1.+1.j, 2.+1.j, 3.+1.j], [ 1.+1.j, 2.+1.j, 3.+1.j, 0.+0.j], [ 2.+1.j, 3.+1.j, 0.+0.j, 0.+0.j], [ 3.+1.j, 0.+0.j, 0.+0.j, 0.+0.j]]) scipy.linalg.toeplitz(np.arange(4)+1.j) Warning: column and row values don't agree; column value used. array([[ 0.+1.j, 1.+1.j, 2.+1.j, 3.+1.j], [ 1.-1.j, 0.+1.j, 1.+1.j, 2.+1.j], [ 2.-1.j, 1.-1.j, 0.+1.j, 1.+1.j], [ 3.-1.j, 2.-1.j, 1.-1.j, 0.+1.j]])
requiring the diagonal to be real might be too tough, I'm not using complex toeplitz, but for other cases, I often have "almost real" values, up to numerical imprecision.
scipy.linalg.toeplitz(np.arange(4)*1.j,-np.arange(4)*1.j) array([[ 0.+0.j, 0.-1.j, 0.-2.j, 0.-3.j], [ 0.+1.j, 0.+0.j, 0.-1.j, 0.-2.j], [ 0.+2.j, 0.+1.j, 0.+0.j, 0.-1.j], [ 0.+3.j, 0.+2.j, 0.+1.j, 0.+0.j]])
also for real values r[0] != c[0] might not be satisfied if there is numerical imprecision and they are only almost equal. Josef
Warren
_______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
josef.pktd@gmail.com wrote:
On Sun, Mar 7, 2010 at 8:32 PM, Warren Weckesser <warren.weckesser@enthought.com> wrote:
Warren Weckesser wrote:
I am going to update linalg.toeplitz to fix the problem reported in ticket #667, but I have some questions about the desired behavior of the function when either `c` (the first column) or `r` (the first row) is complex. The docstring does not say anything about complex values, but the code does some undocumented and unexpected things in this case, one example of which is the bug reported in the ticket #667.
I would like to fix this, but I would first like opinions on exactly how it *should* handle complex arguments. I *think* the idea is that if `c` is complex and `r` is not given, then a Hermitian Toeplitz matrix is created. Is that correct? Note that this requires that c[0] be real, since c[0] is the value that will be in the diagonal of the matrix.
It looks like the behavior was copied from Matlab: http://www.mathworks.com/access/helpdesk/help/techdoc/ref/toeplitz.html
While I am fixing this, I would like to remove the printed warning that occurs when r[0] != c[0]. I would like to change this to raise a ValueError--that is, adopt the philosophy that either the arguments are correct (and the code handles them as [not yet, but soon to be] documented in the docstring, with no warnings needed), or the arguments are wrong and an exception should be raised.
This means that when `r` is not given, an exception would be raised if c[0] is not real.
A very different way to handle the case of `r` not being given is to assume it means a square matrix should be created with r[1:] = 0 (i.e. zeros in all the upper off-diagonals). This avoids the whole issue of complex numbers and conjugates, but it is also a more drastic change.
a lot of code relies on toeplitz of a single argument to build the symmetric matrix.
OK, thanks--that's good to know.
By the way, that is how linalg.hankel handles the case where `r` is not given.
??
I'm not sure what your question marks mean. What I meant by "that is how linalg.hankel handles [it]" was that when `r` is not given, hankel uses zeros to fill in the bottom row. Your examples demonstrate this. My point is moot, since if there is "a lot of code" that relies on creating a symmetric/hermitian matrix, it would be a bad idea to make such a significant change to the API.
import scipy.linalg scipy.linalg.toeplitz(np.arange(4))
array([[0, 1, 2, 3], [1, 0, 1, 2], [2, 1, 0, 1], [3, 2, 1, 0]])
scipy.linalg.hankel(np.arange(4))
array([[ 0., 1., 2., 3.], [ 1., 2., 3., 0.], [ 2., 3., 0., 0.], [ 3., 0., 0., 0.]])
scipy.linalg.hankel(np.arange(4)+1.j)
array([[ 0.+1.j, 1.+1.j, 2.+1.j, 3.+1.j], [ 1.+1.j, 2.+1.j, 3.+1.j, 0.+0.j], [ 2.+1.j, 3.+1.j, 0.+0.j, 0.+0.j], [ 3.+1.j, 0.+0.j, 0.+0.j, 0.+0.j]])
scipy.linalg.toeplitz(np.arange(4)+1.j)
Warning: column and row values don't agree; column value used. array([[ 0.+1.j, 1.+1.j, 2.+1.j, 3.+1.j], [ 1.-1.j, 0.+1.j, 1.+1.j, 2.+1.j], [ 2.-1.j, 1.-1.j, 0.+1.j, 1.+1.j], [ 3.-1.j, 2.-1.j, 1.-1.j, 0.+1.j]])
requiring the diagonal to be real might be too tough, I'm not using complex toeplitz, but for other cases, I often have "almost real" values, up to numerical imprecision.
This could be dealt with like Matlab does: don't worry about it, and just use c[1:] to fill in r[1:]. If c[0] is has nonzero imaginary part, then the matrix is not Hermitian--but that's OK, since that is what the user gave to the function. Your example illustrates one of the bugs: note that, except for the first element, the first column is the *conjugate* of the first argument. The docstring says that `c` is the first column (and the optional second argument is the first row). So it is the elements in the first row starting in the second column that should be conjugated. It also prints that annoying warning.
scipy.linalg.toeplitz(np.arange(4)*1.j,-np.arange(4)*1.j)
array([[ 0.+0.j, 0.-1.j, 0.-2.j, 0.-3.j], [ 0.+1.j, 0.+0.j, 0.-1.j, 0.-2.j], [ 0.+2.j, 0.+1.j, 0.+0.j, 0.-1.j], [ 0.+3.j, 0.+2.j, 0.+1.j, 0.+0.j]])
also for real values r[0] != c[0] might not be satisfied if there is numerical imprecision and they are only almost equal.
True. So if the requirement that r[0] == c[0] is strict, a user who knows the first values are "close enough", but maybe not exactly equal, would have to do something like: r[0] = c[0] t = toeplitz(c, r) I could live with that, but not everyone will agree. Instead, the function could require that the values be "close", using something like numpy.allclose(), but how close is close enough? Should the function then have `atol` and `rtol` arguments to use for the comparison? This feels like overkill, and is not my preference. So that leaves the behavior used by Matlab: use c[0], and ignore r[0]. But I don't like the printed warning (or any logged warning) when c[0] != r[0]. If the program is behaving correctly, and as documented, no warning should be necessary. (If I could start from scratch, I would simply define the `r` argument to give the values from the second column onward, and so completely avoid the issue.) Warren
Two more questions about the toeplitz API: The signature is toeplitz(c, r=None). In the current implementation, if either c or r is a scalar, then c is returned. This results in the following: In [1]: from scipy.linalg import toeplitz In [2]: toeplitz(1, [1,2,3,4]) Out[2]: 1 In [3]: toeplitz([1], [1,2,3,4]) Out[3]: array([[1, 2, 3, 4]]) I would expect these to produce the same result--the second result, in fact. Similarly, the following is surprising: In [4]: toeplitz([1,2,3,4], 1) Out[4]: [1, 2, 3, 4] In [5]: toeplitz([1,2,3,4], [1]) Out[5]: array([[1], [2], [3], [4]]) I think it makes more sense for toeplitz to simply treat a scalar as a 1D array of length 1. Moreover, it seems to me that toeplitz should *always* return a 2D array. My second question: does anyone know why toeplitz uses asarray_chkfinite()? If, in general, numpy arrays can hold the values nan and inf, why should this function care if its values are finite? Any objections to not using asarray_chkfinite()? Warren Warren Weckesser wrote:
josef.pktd@gmail.com wrote:
On Sun, Mar 7, 2010 at 8:32 PM, Warren Weckesser <warren.weckesser@enthought.com> wrote:
Warren Weckesser wrote:
I am going to update linalg.toeplitz to fix the problem reported in ticket #667, but I have some questions about the desired behavior of the function when either `c` (the first column) or `r` (the first row) is complex. The docstring does not say anything about complex values, but the code does some undocumented and unexpected things in this case, one example of which is the bug reported in the ticket #667.
I would like to fix this, but I would first like opinions on exactly how it *should* handle complex arguments. I *think* the idea is that if `c` is complex and `r` is not given, then a Hermitian Toeplitz matrix is created. Is that correct? Note that this requires that c[0] be real, since c[0] is the value that will be in the diagonal of the matrix.
It looks like the behavior was copied from Matlab:
http://www.mathworks.com/access/helpdesk/help/techdoc/ref/toeplitz.html
While I am fixing this, I would like to remove the printed warning that occurs when r[0] != c[0]. I would like to change this to raise a ValueError--that is, adopt the philosophy that either the arguments are correct (and the code handles them as [not yet, but soon to be] documented in the docstring, with no warnings needed), or the arguments are wrong and an exception should be raised.
This means that when `r` is not given, an exception would be raised if c[0] is not real.
A very different way to handle the case of `r` not being given is to assume it means a square matrix should be created with r[1:] = 0 (i.e. zeros in all the upper off-diagonals). This avoids the whole issue of complex numbers and conjugates, but it is also a more drastic change.
a lot of code relies on toeplitz of a single argument to build the symmetric matrix.
OK, thanks--that's good to know.
By the way, that is how linalg.hankel handles the case where `r` is not given.
??
I'm not sure what your question marks mean. What I meant by "that is how linalg.hankel handles [it]" was that when `r` is not given, hankel uses zeros to fill in the bottom row. Your examples demonstrate this. My point is moot, since if there is "a lot of code" that relies on creating a symmetric/hermitian matrix, it would be a bad idea to make such a significant change to the API.
import scipy.linalg scipy.linalg.toeplitz(np.arange(4))
array([[0, 1, 2, 3], [1, 0, 1, 2], [2, 1, 0, 1], [3, 2, 1, 0]])
scipy.linalg.hankel(np.arange(4))
array([[ 0., 1., 2., 3.], [ 1., 2., 3., 0.], [ 2., 3., 0., 0.], [ 3., 0., 0., 0.]])
scipy.linalg.hankel(np.arange(4)+1.j)
array([[ 0.+1.j, 1.+1.j, 2.+1.j, 3.+1.j], [ 1.+1.j, 2.+1.j, 3.+1.j, 0.+0.j], [ 2.+1.j, 3.+1.j, 0.+0.j, 0.+0.j], [ 3.+1.j, 0.+0.j, 0.+0.j, 0.+0.j]])
scipy.linalg.toeplitz(np.arange(4)+1.j)
Warning: column and row values don't agree; column value used. array([[ 0.+1.j, 1.+1.j, 2.+1.j, 3.+1.j], [ 1.-1.j, 0.+1.j, 1.+1.j, 2.+1.j], [ 2.-1.j, 1.-1.j, 0.+1.j, 1.+1.j], [ 3.-1.j, 2.-1.j, 1.-1.j, 0.+1.j]])
requiring the diagonal to be real might be too tough, I'm not using complex toeplitz, but for other cases, I often have "almost real" values, up to numerical imprecision.
This could be dealt with like Matlab does: don't worry about it, and just use c[1:] to fill in r[1:]. If c[0] is has nonzero imaginary part, then the matrix is not Hermitian--but that's OK, since that is what the user gave to the function.
Your example illustrates one of the bugs: note that, except for the first element, the first column is the *conjugate* of the first argument. The docstring says that `c` is the first column (and the optional second argument is the first row). So it is the elements in the first row starting in the second column that should be conjugated. It also prints that annoying warning.
scipy.linalg.toeplitz(np.arange(4)*1.j,-np.arange(4)*1.j)
array([[ 0.+0.j, 0.-1.j, 0.-2.j, 0.-3.j], [ 0.+1.j, 0.+0.j, 0.-1.j, 0.-2.j], [ 0.+2.j, 0.+1.j, 0.+0.j, 0.-1.j], [ 0.+3.j, 0.+2.j, 0.+1.j, 0.+0.j]])
also for real values r[0] != c[0] might not be satisfied if there is numerical imprecision and they are only almost equal.
True. So if the requirement that r[0] == c[0] is strict, a user who knows the first values are "close enough", but maybe not exactly equal, would have to do something like:
r[0] = c[0] t = toeplitz(c, r)
I could live with that, but not everyone will agree.
Instead, the function could require that the values be "close", using something like numpy.allclose(), but how close is close enough? Should the function then have `atol` and `rtol` arguments to use for the comparison? This feels like overkill, and is not my preference.
So that leaves the behavior used by Matlab: use c[0], and ignore r[0]. But I don't like the printed warning (or any logged warning) when c[0] != r[0]. If the program is behaving correctly, and as documented, no warning should be necessary.
(If I could start from scratch, I would simply define the `r` argument to give the values from the second column onward, and so completely avoid the issue.)
Warren
_______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
On Mon, Mar 8, 2010 at 10:25 PM, Warren Weckesser <warren.weckesser@enthought.com> wrote:
Two more questions about the toeplitz API:
I did a (partial) search in trac and on the mailing lists and didn't find any information
The signature is toeplitz(c, r=None). In the current implementation, if either c or r is a scalar, then c is returned. This results in the following:
In [1]: from scipy.linalg import toeplitz
In [2]: toeplitz(1, [1,2,3,4]) Out[2]: 1
In [3]: toeplitz([1], [1,2,3,4]) Out[3]: array([[1, 2, 3, 4]])
I would expect these to produce the same result--the second result, in fact.
Similarly, the following is surprising:
In [4]: toeplitz([1,2,3,4], 1) Out[4]: [1, 2, 3, 4]
In [5]: toeplitz([1,2,3,4], [1]) Out[5]: array([[1], [2], [3], [4]])
I think it makes more sense for toeplitz to simply treat a scalar as a 1D array of length 1. Moreover, it seems to me that toeplitz should *always* return a 2D array.
np.kron(5,3) 15 np.kron([5],[3]) array([15]) np.kron([5],[3]).shape (1,) scipy.linalg.hankel(1) 1 scipy.linalg.hankel(1,2) 1 np.triu(3) Traceback (most recent call last): File "C:\Programs\Python25\Lib\site-packages\numpy\lib\twodim_base.py",
np.triu([3]) Traceback (most recent call last): File "C:\Programs\Python25\Lib\site-packages\numpy\lib\twodim_base.py",
numpy/scipy is not matlab, so the policy is still a bit vague line 442, in triu out = multiply((1-tri(m.shape[0], m.shape[1], k-1, int)),m) IndexError: tuple index out of range line 442, in triu out = multiply((1-tri(m.shape[0], m.shape[1], k-1, int)),m) IndexError: tuple index out of range
np.triu([[3]]) array([[3]])
scipy.linalg.kron([-1j], [1.j]) Traceback (most recent call last): File "\Programs\Python25\Lib\site-packages\scipy\linalg\basic.py", line 835, in kron AttributeError: 'list' object has no attribute 'flags'
and I filed a ticket for scipy.linalg.block_diag because there I can see a use for scalar to 2d. I don't have a strong opinion about always 2d result, but a bit more consistency in numpy/scipy would be useful.
My second question: does anyone know why toeplitz uses asarray_chkfinite()? If, in general, numpy arrays can hold the values nan and inf, why should this function care if its values are finite? Any objections to not using asarray_chkfinite()?
I agree, I don't think it's the job of toeplitz or hankle to throw an exception if I want to have a inf or nan in my matrix.
Warren
Warren Weckesser wrote:
josef.pktd@gmail.com wrote:
On Sun, Mar 7, 2010 at 8:32 PM, Warren Weckesser <warren.weckesser@enthought.com> wrote:
Warren Weckesser wrote:
I am going to update linalg.toeplitz to fix the problem reported in ticket #667, but I have some questions about the desired behavior of the function when either `c` (the first column) or `r` (the first row) is complex. The docstring does not say anything about complex values, but the code does some undocumented and unexpected things in this case, one example of which is the bug reported in the ticket #667.
I would like to fix this, but I would first like opinions on exactly how it *should* handle complex arguments. I *think* the idea is that if `c` is complex and `r` is not given, then a Hermitian Toeplitz matrix is created. Is that correct? Note that this requires that c[0] be real, since c[0] is the value that will be in the diagonal of the matrix.
It looks like the behavior was copied from Matlab:
http://www.mathworks.com/access/helpdesk/help/techdoc/ref/toeplitz.html
While I am fixing this, I would like to remove the printed warning that occurs when r[0] != c[0]. I would like to change this to raise a ValueError--that is, adopt the philosophy that either the arguments are correct (and the code handles them as [not yet, but soon to be] documented in the docstring, with no warnings needed), or the arguments are wrong and an exception should be raised.
This means that when `r` is not given, an exception would be raised if c[0] is not real.
A very different way to handle the case of `r` not being given is to assume it means a square matrix should be created with r[1:] = 0 (i.e. zeros in all the upper off-diagonals). This avoids the whole issue of complex numbers and conjugates, but it is also a more drastic change.
a lot of code relies on toeplitz of a single argument to build the symmetric matrix.
OK, thanks--that's good to know.
By the way, that is how linalg.hankel handles the case where `r` is not given.
??
I'm not sure what your question marks mean. What I meant by "that is how linalg.hankel handles [it]" was that when `r` is not given, hankel uses zeros to fill in the bottom row. Your examples demonstrate this. My point is moot, since if there is "a lot of code" that relies on creating a symmetric/hermitian matrix, it would be a bad idea to make such a significant change to the API.
import scipy.linalg scipy.linalg.toeplitz(np.arange(4))
array([[0, 1, 2, 3], [1, 0, 1, 2], [2, 1, 0, 1], [3, 2, 1, 0]])
scipy.linalg.hankel(np.arange(4))
array([[ 0., 1., 2., 3.], [ 1., 2., 3., 0.], [ 2., 3., 0., 0.], [ 3., 0., 0., 0.]])
scipy.linalg.hankel(np.arange(4)+1.j)
array([[ 0.+1.j, 1.+1.j, 2.+1.j, 3.+1.j], [ 1.+1.j, 2.+1.j, 3.+1.j, 0.+0.j], [ 2.+1.j, 3.+1.j, 0.+0.j, 0.+0.j], [ 3.+1.j, 0.+0.j, 0.+0.j, 0.+0.j]])
scipy.linalg.toeplitz(np.arange(4)+1.j)
Warning: column and row values don't agree; column value used. array([[ 0.+1.j, 1.+1.j, 2.+1.j, 3.+1.j], [ 1.-1.j, 0.+1.j, 1.+1.j, 2.+1.j], [ 2.-1.j, 1.-1.j, 0.+1.j, 1.+1.j], [ 3.-1.j, 2.-1.j, 1.-1.j, 0.+1.j]])
requiring the diagonal to be real might be too tough, I'm not using complex toeplitz, but for other cases, I often have "almost real" values, up to numerical imprecision.
This could be dealt with like Matlab does: don't worry about it, and just use c[1:] to fill in r[1:]. If c[0] is has nonzero imaginary part, then the matrix is not Hermitian--but that's OK, since that is what the user gave to the function.
Your example illustrates one of the bugs: note that, except for the first element, the first column is the *conjugate* of the first argument. The docstring says that `c` is the first column (and the optional second argument is the first row). So it is the elements in the first row starting in the second column that should be conjugated. It also prints that annoying warning.
That's how I interpreted the bug
scipy.linalg.toeplitz(np.arange(4)*1.j,-np.arange(4)*1.j)
array([[ 0.+0.j, 0.-1.j, 0.-2.j, 0.-3.j], [ 0.+1.j, 0.+0.j, 0.-1.j, 0.-2.j], [ 0.+2.j, 0.+1.j, 0.+0.j, 0.-1.j], [ 0.+3.j, 0.+2.j, 0.+1.j, 0.+0.j]])
also for real values r[0] != c[0] might not be satisfied if there is numerical imprecision and they are only almost equal.
True. So if the requirement that r[0] == c[0] is strict, a user who knows the first values are "close enough", but maybe not exactly equal, would have to do something like:
r[0] = c[0] t = toeplitz(c, r)
I could live with that, but not everyone will agree.
Instead, the function could require that the values be "close", using something like numpy.allclose(), but how close is close enough? Should the function then have `atol` and `rtol` arguments to use for the comparison? This feels like overkill, and is not my preference.
So that leaves the behavior used by Matlab: use c[0], and ignore r[0]. But I don't like the printed warning (or any logged warning) when c[0] != r[0]. If the program is behaving correctly, and as documented, no warning should be necessary.
raising a proper Warning might be useful and avoid the almost equal issue I think, all the uses I have seen use only a single argument to build the symmetric toeplitz matrix. Josef
(If I could start from scratch, I would simply define the `r` argument to give the values from the second column onward, and so completely avoid the issue.)
Warren
_______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
_______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
josef.pktd@gmail.com wrote:
On Mon, Mar 8, 2010 at 10:25 PM, Warren Weckesser <warren.weckesser@enthought.com> wrote:
Two more questions about the toeplitz API:
I did a (partial) search in trac and on the mailing lists and didn't find any information
I have submitted a patch, attached to ticket #667. The patch changes the behavior of toeplitz() and hankel(), and adds a new function, circulant(). The patch also adds several tests of these functions. I can commit the patch, but I'd like to know if y'all agree with the changes first. Warren
The signature is toeplitz(c, r=None). In the current implementation, if either c or r is a scalar, then c is returned. This results in the following:
In [1]: from scipy.linalg import toeplitz
In [2]: toeplitz(1, [1,2,3,4]) Out[2]: 1
In [3]: toeplitz([1], [1,2,3,4]) Out[3]: array([[1, 2, 3, 4]])
I would expect these to produce the same result--the second result, in fact.
Similarly, the following is surprising:
In [4]: toeplitz([1,2,3,4], 1) Out[4]: [1, 2, 3, 4]
In [5]: toeplitz([1,2,3,4], [1]) Out[5]: array([[1], [2], [3], [4]])
I think it makes more sense for toeplitz to simply treat a scalar as a 1D array of length 1. Moreover, it seems to me that toeplitz should *always* return a 2D array.
numpy/scipy is not matlab, so the policy is still a bit vague
np.kron(5,3)
15
np.kron([5],[3])
array([15])
np.kron([5],[3]).shape
(1,)
scipy.linalg.hankel(1)
1
scipy.linalg.hankel(1,2)
1
np.triu(3)
Traceback (most recent call last): File "C:\Programs\Python25\Lib\site-packages\numpy\lib\twodim_base.py", line 442, in triu out = multiply((1-tri(m.shape[0], m.shape[1], k-1, int)),m) IndexError: tuple index out of range
np.triu([3])
Traceback (most recent call last): File "C:\Programs\Python25\Lib\site-packages\numpy\lib\twodim_base.py", line 442, in triu out = multiply((1-tri(m.shape[0], m.shape[1], k-1, int)),m) IndexError: tuple index out of range
np.triu([[3]])
array([[3]])
scipy.linalg.kron([-1j], [1.j])
Traceback (most recent call last): File "\Programs\Python25\Lib\site-packages\scipy\linalg\basic.py", line 835, in kron AttributeError: 'list' object has no attribute 'flags'
and I filed a ticket for scipy.linalg.block_diag because there I can see a use for scalar to 2d.
I don't have a strong opinion about always 2d result, but a bit more consistency in numpy/scipy would be useful.
My second question: does anyone know why toeplitz uses asarray_chkfinite()? If, in general, numpy arrays can hold the values nan and inf, why should this function care if its values are finite? Any objections to not using asarray_chkfinite()?
I agree, I don't think it's the job of toeplitz or hankle to throw an exception if I want to have a inf or nan in my matrix.
Warren
Warren Weckesser wrote:
josef.pktd@gmail.com wrote:
On Sun, Mar 7, 2010 at 8:32 PM, Warren Weckesser <warren.weckesser@enthought.com> wrote:
Warren Weckesser wrote:
I am going to update linalg.toeplitz to fix the problem reported in ticket #667, but I have some questions about the desired behavior of the function when either `c` (the first column) or `r` (the first row) is complex. The docstring does not say anything about complex values, but the code does some undocumented and unexpected things in this case, one example of which is the bug reported in the ticket #667.
I would like to fix this, but I would first like opinions on exactly how it *should* handle complex arguments. I *think* the idea is that if `c` is complex and `r` is not given, then a Hermitian Toeplitz matrix is created. Is that correct? Note that this requires that c[0] be real, since c[0] is the value that will be in the diagonal of the matrix.
It looks like the behavior was copied from Matlab:
http://www.mathworks.com/access/helpdesk/help/techdoc/ref/toeplitz.html
While I am fixing this, I would like to remove the printed warning that occurs when r[0] != c[0]. I would like to change this to raise a ValueError--that is, adopt the philosophy that either the arguments are correct (and the code handles them as [not yet, but soon to be] documented in the docstring, with no warnings needed), or the arguments are wrong and an exception should be raised.
This means that when `r` is not given, an exception would be raised if c[0] is not real.
A very different way to handle the case of `r` not being given is to assume it means a square matrix should be created with r[1:] = 0 (i.e. zeros in all the upper off-diagonals). This avoids the whole issue of complex numbers and conjugates, but it is also a more drastic change.
a lot of code relies on toeplitz of a single argument to build the symmetric matrix.
OK, thanks--that's good to know.
By the way, that is how linalg.hankel handles the case where `r` is not given.
??
I'm not sure what your question marks mean. What I meant by "that is how linalg.hankel handles [it]" was that when `r` is not given, hankel uses zeros to fill in the bottom row. Your examples demonstrate this. My point is moot, since if there is "a lot of code" that relies on creating a symmetric/hermitian matrix, it would be a bad idea to make such a significant change to the API.
> import scipy.linalg > scipy.linalg.toeplitz(np.arange(4)) > > > array([[0, 1, 2, 3], [1, 0, 1, 2], [2, 1, 0, 1], [3, 2, 1, 0]])
> scipy.linalg.hankel(np.arange(4)) > > > array([[ 0., 1., 2., 3.], [ 1., 2., 3., 0.], [ 2., 3., 0., 0.], [ 3., 0., 0., 0.]])
> scipy.linalg.hankel(np.arange(4)+1.j) > > > array([[ 0.+1.j, 1.+1.j, 2.+1.j, 3.+1.j], [ 1.+1.j, 2.+1.j, 3.+1.j, 0.+0.j], [ 2.+1.j, 3.+1.j, 0.+0.j, 0.+0.j], [ 3.+1.j, 0.+0.j, 0.+0.j, 0.+0.j]])
> scipy.linalg.toeplitz(np.arange(4)+1.j) > > > Warning: column and row values don't agree; column value used. array([[ 0.+1.j, 1.+1.j, 2.+1.j, 3.+1.j], [ 1.-1.j, 0.+1.j, 1.+1.j, 2.+1.j], [ 2.-1.j, 1.-1.j, 0.+1.j, 1.+1.j], [ 3.-1.j, 2.-1.j, 1.-1.j, 0.+1.j]])
requiring the diagonal to be real might be too tough, I'm not using complex toeplitz, but for other cases, I often have "almost real" values, up to numerical imprecision.
This could be dealt with like Matlab does: don't worry about it, and just use c[1:] to fill in r[1:]. If c[0] is has nonzero imaginary part, then the matrix is not Hermitian--but that's OK, since that is what the user gave to the function.
Your example illustrates one of the bugs: note that, except for the first element, the first column is the *conjugate* of the first argument. The docstring says that `c` is the first column (and the optional second argument is the first row). So it is the elements in the first row starting in the second column that should be conjugated. It also prints that annoying warning.
That's how I interpreted the bug
> scipy.linalg.toeplitz(np.arange(4)*1.j,-np.arange(4)*1.j) > > > array([[ 0.+0.j, 0.-1.j, 0.-2.j, 0.-3.j], [ 0.+1.j, 0.+0.j, 0.-1.j, 0.-2.j], [ 0.+2.j, 0.+1.j, 0.+0.j, 0.-1.j], [ 0.+3.j, 0.+2.j, 0.+1.j, 0.+0.j]])
also for real values r[0] != c[0] might not be satisfied if there is numerical imprecision and they are only almost equal.
True. So if the requirement that r[0] == c[0] is strict, a user who knows the first values are "close enough", but maybe not exactly equal, would have to do something like:
r[0] = c[0] t = toeplitz(c, r)
I could live with that, but not everyone will agree.
Instead, the function could require that the values be "close", using something like numpy.allclose(), but how close is close enough? Should the function then have `atol` and `rtol` arguments to use for the comparison? This feels like overkill, and is not my preference.
So that leaves the behavior used by Matlab: use c[0], and ignore r[0]. But I don't like the printed warning (or any logged warning) when c[0] != r[0]. If the program is behaving correctly, and as documented, no warning should be necessary.
raising a proper Warning might be useful and avoid the almost equal issue
I think, all the uses I have seen use only a single argument to build the symmetric toeplitz matrix.
Josef
(If I could start from scratch, I would simply define the `r` argument to give the values from the second column onward, and so completely avoid the issue.)
Warren
_______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
_______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
_______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
On Sat, Mar 20, 2010 at 9:29 PM, Warren Weckesser <warren.weckesser@enthought.com> wrote:
josef.pktd@gmail.com wrote:
On Mon, Mar 8, 2010 at 10:25 PM, Warren Weckesser <warren.weckesser@enthought.com> wrote:
Two more questions about the toeplitz API:
I did a (partial) search in trac and on the mailing lists and didn't find any information
I have submitted a patch, attached to ticket #667. The patch changes the behavior of toeplitz() and hankel(), and adds a new function, circulant(). The patch also adds several tests of these functions. I can commit the patch, but I'd like to know if y'all agree with the changes first.
fine with me Josef
Warren
The signature is toeplitz(c, r=None). In the current implementation, if either c or r is a scalar, then c is returned. This results in the following:
In [1]: from scipy.linalg import toeplitz
In [2]: toeplitz(1, [1,2,3,4]) Out[2]: 1
In [3]: toeplitz([1], [1,2,3,4]) Out[3]: array([[1, 2, 3, 4]])
I would expect these to produce the same result--the second result, in fact.
Similarly, the following is surprising:
In [4]: toeplitz([1,2,3,4], 1) Out[4]: [1, 2, 3, 4]
In [5]: toeplitz([1,2,3,4], [1]) Out[5]: array([[1], [2], [3], [4]])
I think it makes more sense for toeplitz to simply treat a scalar as a 1D array of length 1. Moreover, it seems to me that toeplitz should *always* return a 2D array.
numpy/scipy is not matlab, so the policy is still a bit vague
np.kron(5,3)
15
np.kron([5],[3])
array([15])
np.kron([5],[3]).shape
(1,)
scipy.linalg.hankel(1)
1
scipy.linalg.hankel(1,2)
1
np.triu(3)
Traceback (most recent call last): File "C:\Programs\Python25\Lib\site-packages\numpy\lib\twodim_base.py", line 442, in triu out = multiply((1-tri(m.shape[0], m.shape[1], k-1, int)),m) IndexError: tuple index out of range
np.triu([3])
Traceback (most recent call last): File "C:\Programs\Python25\Lib\site-packages\numpy\lib\twodim_base.py", line 442, in triu out = multiply((1-tri(m.shape[0], m.shape[1], k-1, int)),m) IndexError: tuple index out of range
np.triu([[3]])
array([[3]])
scipy.linalg.kron([-1j], [1.j])
Traceback (most recent call last): File "\Programs\Python25\Lib\site-packages\scipy\linalg\basic.py", line 835, in kron AttributeError: 'list' object has no attribute 'flags'
and I filed a ticket for scipy.linalg.block_diag because there I can see a use for scalar to 2d.
I don't have a strong opinion about always 2d result, but a bit more consistency in numpy/scipy would be useful.
My second question: does anyone know why toeplitz uses asarray_chkfinite()? If, in general, numpy arrays can hold the values nan and inf, why should this function care if its values are finite? Any objections to not using asarray_chkfinite()?
I agree, I don't think it's the job of toeplitz or hankle to throw an exception if I want to have a inf or nan in my matrix.
Warren
Warren Weckesser wrote:
josef.pktd@gmail.com wrote:
On Sun, Mar 7, 2010 at 8:32 PM, Warren Weckesser <warren.weckesser@enthought.com> wrote:
Warren Weckesser wrote:
> I am going to update linalg.toeplitz to fix the problem reported in > ticket #667, but I have some questions about the desired behavior of > the function when either `c` (the first column) or `r` (the first row) > is complex. The docstring does not say anything about complex values, > but the code does some undocumented and unexpected things in this > case, one example of which is the bug reported in the ticket #667. > > I would like to fix this, but I would first like opinions on exactly > how it *should* handle complex arguments. I *think* the idea is that > if `c` is complex and `r` is not given, then a Hermitian Toeplitz > matrix is created. Is that correct? Note that this requires that > c[0] be real, since c[0] is the value that will be in the diagonal of > the matrix. > > > > It looks like the behavior was copied from Matlab:
http://www.mathworks.com/access/helpdesk/help/techdoc/ref/toeplitz.html
> While I am fixing this, I would like to remove the printed warning > that occurs when r[0] != c[0]. I would like to change this to raise > a ValueError--that is, adopt the philosophy that either the arguments > are correct (and the code handles them as [not yet, but soon to be] > documented in the docstring, with no warnings needed), or the > arguments are wrong and an exception should be raised. > > This means that when `r` is not given, an exception would be raised > if c[0] is not real. > > A very different way to handle the case of `r` not being given is to > assume it means a square matrix should be created with r[1:] = 0 > (i.e. zeros in all the upper off-diagonals). This avoids the whole > issue of complex numbers and conjugates, but it is also a more > drastic change. > > > a lot of code relies on toeplitz of a single argument to build the symmetric matrix.
OK, thanks--that's good to know.
> By the way, that is how linalg.hankel handles the case where `r` is not given.
??
I'm not sure what your question marks mean. What I meant by "that is how linalg.hankel handles [it]" was that when `r` is not given, hankel uses zeros to fill in the bottom row. Your examples demonstrate this. My point is moot, since if there is "a lot of code" that relies on creating a symmetric/hermitian matrix, it would be a bad idea to make such a significant change to the API.
>> import scipy.linalg >> scipy.linalg.toeplitz(np.arange(4)) >> >> >> array([[0, 1, 2, 3], [1, 0, 1, 2], [2, 1, 0, 1], [3, 2, 1, 0]])
>> scipy.linalg.hankel(np.arange(4)) >> >> >> array([[ 0., 1., 2., 3.], [ 1., 2., 3., 0.], [ 2., 3., 0., 0.], [ 3., 0., 0., 0.]])
>> scipy.linalg.hankel(np.arange(4)+1.j) >> >> >> array([[ 0.+1.j, 1.+1.j, 2.+1.j, 3.+1.j], [ 1.+1.j, 2.+1.j, 3.+1.j, 0.+0.j], [ 2.+1.j, 3.+1.j, 0.+0.j, 0.+0.j], [ 3.+1.j, 0.+0.j, 0.+0.j, 0.+0.j]])
>> scipy.linalg.toeplitz(np.arange(4)+1.j) >> >> >> Warning: column and row values don't agree; column value used. array([[ 0.+1.j, 1.+1.j, 2.+1.j, 3.+1.j], [ 1.-1.j, 0.+1.j, 1.+1.j, 2.+1.j], [ 2.-1.j, 1.-1.j, 0.+1.j, 1.+1.j], [ 3.-1.j, 2.-1.j, 1.-1.j, 0.+1.j]])
requiring the diagonal to be real might be too tough, I'm not using complex toeplitz, but for other cases, I often have "almost real" values, up to numerical imprecision.
This could be dealt with like Matlab does: don't worry about it, and just use c[1:] to fill in r[1:]. If c[0] is has nonzero imaginary part, then the matrix is not Hermitian--but that's OK, since that is what the user gave to the function.
Your example illustrates one of the bugs: note that, except for the first element, the first column is the *conjugate* of the first argument. The docstring says that `c` is the first column (and the optional second argument is the first row). So it is the elements in the first row starting in the second column that should be conjugated. It also prints that annoying warning.
That's how I interpreted the bug
>> scipy.linalg.toeplitz(np.arange(4)*1.j,-np.arange(4)*1.j) >> >> >> array([[ 0.+0.j, 0.-1.j, 0.-2.j, 0.-3.j], [ 0.+1.j, 0.+0.j, 0.-1.j, 0.-2.j], [ 0.+2.j, 0.+1.j, 0.+0.j, 0.-1.j], [ 0.+3.j, 0.+2.j, 0.+1.j, 0.+0.j]])
also for real values r[0] != c[0] might not be satisfied if there is numerical imprecision and they are only almost equal.
True. So if the requirement that r[0] == c[0] is strict, a user who knows the first values are "close enough", but maybe not exactly equal, would have to do something like:
r[0] = c[0] t = toeplitz(c, r)
I could live with that, but not everyone will agree.
Instead, the function could require that the values be "close", using something like numpy.allclose(), but how close is close enough? Should the function then have `atol` and `rtol` arguments to use for the comparison? This feels like overkill, and is not my preference.
So that leaves the behavior used by Matlab: use c[0], and ignore r[0]. But I don't like the printed warning (or any logged warning) when c[0] != r[0]. If the program is behaving correctly, and as documented, no warning should be necessary.
raising a proper Warning might be useful and avoid the almost equal issue
I think, all the uses I have seen use only a single argument to build the symmetric toeplitz matrix.
Josef
(If I could start from scratch, I would simply define the `r` argument to give the values from the second column onward, and so completely avoid the issue.)
Warren
_______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
_______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
_______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
_______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
participants (2)
-
josef.pktd＠gmail.com -
Warren Weckesser