Arithmetic operations between QT matrices
Contents
Additions and multiplications
If we define two QT matrices $A$ and $B$, we can easily compute $C = A \pm B$ using the familiar matrix operations in MATLAB.
A = cqt([ 4 2 1 ], [4 -1 2]); B = cqt([ 6 3 1 1 ], [ 6 1 ]); C = A + B
C = CQT Matrix of size Inf x Inf - Toeplitz part (leading 6 x 5 block): 10 0 2 0 0 5 10 0 2 0 2 5 10 0 2 1 2 5 10 0 0 1 2 5 10 0 0 1 2 5
Matrix multiplication can be handled in a similar way. Note however that in this case the bandwidth of the result is larger, since the symbol $c(z)$ of the product $C = AB$ is the product of the symbols of $A$ and $B$, and therefore has higher degree.
C = A * B
C = CQT Matrix of size Inf x Inf Rank of top-left correction: 1 - Toeplitz part (leading 8 x 6 block): 25.0000 4.0000 11.0000 2.0000 0 0 26.0000 25.0000 4.0000 11.0000 2.0000 0 15.0000 26.0000 25.0000 4.0000 11.0000 2.0000 9.0000 15.0000 26.0000 25.0000 4.0000 11.0000 3.0000 9.0000 15.0000 26.0000 25.0000 4.0000 1.0000 3.0000 9.0000 15.0000 26.0000 25.0000 0 1.0000 3.0000 9.0000 15.0000 26.0000 0 0 1.0000 3.0000 9.0000 15.0000 - Finite correction (top-left corner): -2.0000 -1.0000
Computing inverses
The inv command is overloaded. Although it is rarely a good idea to compute inverses directly, it may sometimes be useful. If a QT matrix is invertible, then its inverse is again in the QT class. The invertibility is equivalent (for infinite matrices) to asking that $a(z) \neq 0$ for $|z| = 1$. This condition is satisfied by both $A$ and $B$ defined above, so we may compute their inverses.
iA = inv(A);
It's importante to notice that in this case the symbol of iA$ is not of finite length, since it is the Laurent expansion of $a(z)^{-1}$. In practice, the toolbox automatically truncate the symbol to the necessary length. We can visualize the exponential decay of the symbol entries by plotting a finite section of the inverse.
surf(abs(iA(1:30, 1:30))); set(gca,'zscale','log'); % log-scale to appreciate the exponential decay title(gca, 'Decay in the entries of A^{-1}');
To check how long is the symbol, or how large is the support of the low-rank correction, we may use the command symbol and correction:
[am, ap] = symbol(iA); nm = length(am), np = length(ap), nc = size(correction(iA))
nm = 35 np = 63 nc = 34 62
These lengths are not symmetric because the symbol of $A$ was not symmetric in the first place.
Linear systems
Quite often, we do not need to explicitly compute the inverse, but rather to solve a linear system of the form
$$ Ax = b, $$
where $A$ is QT, and $x$ and $b$ vectors of infinite length. If this is the case, we may resort to computing a $UL$ factorization instead, and solve the system by running the following commands.
b = randn(4, 1); b = cqt(b); % We interpret b as a vector of infite size. x = A \ b; x = correction(x) % Automatically truncate x to the necessary length
x = 0.1268 -0.0019 0.1081 -0.1005 0.0295 0.0061 -0.0096 0.0037 0.0001 -0.0009 0.0004 -0.0000 -0.0001 0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000 0.0000 0.0000 -0.0000 0.0000 0.0000 -0.0000 0.0000 0.0000 -0.0000 0.0000 0.0000 -0.0000 0.0000
Please visit the chapter on linear systems for further details on this matter.