rluce / tlcomp Goto Github PK

The current ones can be interpreted, but could be made nicer

>> TM1

TM1 = 

  8×8 ToepMat array with properties:

    c: [8×1 double]
    r: [8×1 double]

>> TM1^2

ans = 

  8×8 TLMat array with properties:

    G: [8×4 double]
    B: [8×4 double]

Currently only square matrices are supported in |ToepMat| and |TLMat|. This restriction could be lifted; it will be a lot of work though.

The function |toeplksolvetoeplk| currently computes inv(TL1) * TL2, and this is used in both |mrdivide| and |mldivide| in |ToepMat|. However, the current implementation needs to transpose a lot, e.g.,

TL1 / TL2 = (TL2' \ TL1' )'

which makes generators unneccessarily long, and introduces additional roundoff from the explicit ctransposes involved.

Instead |toeplksolvetoeplk| could be extended so that both operands are used implicitly as ctransposed quantities. This change entails that the matrices involved need to be represented as Zp/Zm or Zm/Zp low-drank, which requires some carful thoughts upfront. One can draw, however, from the already existing functionality in |toeplksolve| and |toeplkmult|.

Polynomials of Topelitz matrices are currently computed by first gathering all necessary mononials, and then adding them up. A Horner scheme should have the same advantages as in the dense, unstructured matrix case.

This should be tested/confirmed, and then offered as an alternative in toeppolyvalm.m.

For |TLMat| and |ToepMat| the matrix 1/inf/fro norms can simply computed in O(n^2) by

norm(full(T), p).

Note that this is done in O(n^2) space. This can easily be brought down to O(n) space and O(n) time for |ToepMat|. For |TLMat|, these operations can be implemented in O(n) space as follows. The columns of T are generated one-after-one in toeplkreconstruct, and once they are accounted for in the norm computation, they can be discarded.

When transposing a TL matrix, the rank increases if the first column / last row is not in some range of the generator. It could well be that this can be checked by normalizing the generator (aka proper form) by hyperbolic rotations.

Moreover, generating the first column is an important operation by itself; it is used for example for norm computation, or, more generally, in any column-traversal based operation.

We want to compute det(TLMat). This can be done by computing the LU factorization after transformation to Cauchy like, and to multiply all the determinants of the unitiary and trianuglar matrices involved. Of course none of the matrices should be computed explicitly, so that the memory footprint is O(n).

Assume given a dense Toeplitz matrix A, then it would be nice to convert it to a |ToepMat| object by

ToepMat(A)

instead the cumbersome ToepMat(A(:,1), A(1,:)).

Even if A is not Toeplitz, the constructor could return the closest Toeplitz matrix to A in the fro norm. Note that this optimization problem is solved by taking means along the diagnoals of A.

To generate a TLMat from a given ToepMat T, one has to write

TLMat(T.c, T.r)

but it would be much nicer to simply write

TLMat(T).

For experimentation it would be nice to have an experimental connection to a superfast HSS solver.

We could easily reduce the dependency on drsolve to a minimum by using "just" the main function clsolve. The basic FFT based pre- and post transformations could be done by us. Some funcionality of these (nroots1, ttimes) is already implemented in other forms and can be refactored.

• ftimes.m
• nroots1.m
• t2cl.m
• tl2cl.m
• tlsolve.m
• tsolve.m
• ttimes.m

For a square, singular matrix A matlab issues warnings if A\b (or b/A) is tried. The behaviour of such singular cases is currently not consistently implemented in \ or / in |TLMat| or |ToepMat|.

Right now we store G, B in TLMat, but storing orthogonal bases separately would allow for more convenient operations.

What should

TM(:,1:5) = A

do?

Some algorithmic choices should actually be user controllable, for example the number of iterative refinement steps, or the kind of multiplication should use full or fft mode. Looking ahead, a lot more options can be expected:

• truncation tolerance
• solver choice (fast/superfast)
• regularization eps for LDU

The option mechanism should supprt both global and local influence. Nice syntax ideas:

global_opt = TLMat.getopts % Global option set, modifiable from outside
TLMat.opt.num_iterref = 1; % Global setting for all TLMat objects
T = tleye(8);
local_opt = T.getopts % Option set pertaining to T
T.opt.num_iterref = 0; % Local setting overrides global one

For a TLMat objectTL , length(TL) should behave just like length(full(TL)), but currently it just returns 1, which doesn't make sense.

Currently istoeplitz checks whether the operand is a Toeplitz matrix by constructing a full toeplitz matrix using toeplitz, and comparing the two. For large matrices it may be faster to iterate over the columns or row of the operand, comparing the cyclically shifted values directly. At minimum, it avoids another O(n^2) memory load.

At every recompression we obtain the |GB'|_2 for free. Since |A|_2 <= sqrt(n) * |GB'|_2, this gives a good norm estimate at zero computational overhead.

Given a matrix A that is potentially close to a TL matrix, we would like to generate a TL matrix by

T = TLMat(A)

so that an efficient representation of A is obtained. This is currently achieved by computing its displacement D(A), and computing a SVD of D(A), which is expensive (cubic operation).

A more economic choice would be to use Adaptive Cross Approximation on D(A) to recover a low rank approximation (essentially d * n^2, where d is the rank). We then should also support syntax for controlling the rank, i.e.

T = TLMat(A, d) % prescribed drank d
T = TLMat(A, tol) % prescribed tolerance pivot_k / pivot_1 that indicates to stop

It is trivial to detect wether a Toeplitz matrix is triangular. We should take advantage of this when

• T\b
• det(T)
• T*b

As an aside we should implement the functions istril and istriu and use them to dispatch the operations.

The displacement rank of a TL matrix A is the rank of its displacement D(A). Currently, whenever computations with the generator D(A) = G*B' are carried out, we compress G,B back to its numerical rank, essentially determined by Matlab's default rank tolerance.

It can be shown that the truncation error in the generator is amplified only by |A|, so that when the generator is truncated to some larger tolerance, the matrix it represents still is close to A. Thus it would be nice to give the user control of the truncation tolerance, say, in the constructor

T = TLMat(G, B, 'eps_trunc', tol)

or even letting her or him control it on the fly

T.set_eps_trunc(new_tol).

The latter makes a lot of sense in an iterative approximation scheme, where the initial accuracy of the approximants are quite bad anyway.

At current state both TLMat and ToepMat are sitting in one monolithic file. Matlab's indention rules (class/methods/fuction) induce a lot of whitespace noise in the code, making it difficult to maintain the 80 char line width limit.

Instead the two classes should sit in the own directories, with each method becoming its own file. This is described in the Matlab documentation for OO development.

There is no least squares capabilities right now. Choice of algorithms is not trivial.

A useful addtion would be a chol function for ToepMat. This would be based on the Schur algorithm. In addition, a cholsolve method should be implemented, that uses the Schur code with implicit triangular solve, resulting in a O(n) memory solver.

TM/TM is currently unsupported. In contrast TM\TM is supported so we could reduce it to the latter. In would involved, however, in a ctranspose(TL), which could possibly be avoided.

T^s is currently not implemented although the corresponding toeppowers is ready to go.

For TL this could/should be implemented only on the TL*TL syntax, unless a compelling reason compes up to defer the computation on the generator leve.