## Using Methods on Operators- Multiplication
- Transposition and conjugation
- Addition and subtraction
- Operator Information
- Dictionaries and arrays
- Block diagonal operators
- Kronecker products
- Subset assignment and reference
- Elementwise operations
- Power, inverse, and backslash operations
- Solving systems of linear equations
- Application to non-linear operators
All Spot operators provide routines for multiplying vectors by itself and its adjoint. The real power of the operator toolbox is the ease with which basic operators (see the Index of Operators) can be combined into more powerful operators, which themselves can be combined into yet more powerful operators. In this section we discuss the type of operations can be applied to all Spot operators. All operations discussed in this section are implemented using underlying meta-operators. Pointers to these operators are given at the end of each section, where relevant. In the cases where the meta-operators provide more extensive possibilities, they are described in more detail below the basic use. We denote by – arbitrary Spot operators and by any matrix. All operator and matrix are assumed to have appropriate dimensions that can vary line per line. Most operations described in this section are meaningful only when applied to linear operators; see Application to non-linear operators for more details. ## MultiplicationMultiplication by operators is the most elementary operation in Spot. In its simplest form we have operator-vector products such as: y = A * x; z = B * y; n = m'* C; This will evaluate the application of on , and then to (the third command is evaluated as ). In both cases the result will be an explicit vector of real or complex numbers. An equivalent approach is to first create a new operator representing the product , and then apply this to to immediately get vector : C = B * A; % Construct compound operator z = C * x; % Evaluate z = BAx Although Spot operators are expected to support only multiplication of vectors, it is possible to write operator-matrix multiplications. However, it should be kept in mind that this is implemented as a series of operator-vector multiplications. The second command below is evaluated as . N = C * M; N = M * C; % A special case of multiplication is multiplication by scalars. Generally, this will give a new operator scaled by that quantity; C = 3 * A; % Construct compound operator C = 3A C = A * 3; One notable exception is when the corresponding dimension of the
operator is one. In that case we have a valid matrix-vector product
which results in a numeric solution. In other words, matrix-vector
products take precedence over scalar multiplication. In such cases it
is still possible to scale the operator, either by changing the order
(as shown above) or by directly setting the C = -A; B = +A; This works regardless of the dimensions of the operator. Note that
elementwise multiplication using In the evaluation and comparison of algorithms it is often useful to know how many multiplications have been performed using a given operator. To facilitate these kind of measurements it is possible to associate a counter variable to each operator.
## Transposition and conjugationAs mentioned in the introduction, each Spot operator implements
multiplication by itself and its conjugate. Using Matlab’s transpose
operator B = A'; x = A'* y; x = B * y; % Identical result as previous line When using transposes within calculations these new operators are
discarded as soon as the multiplication has been done. Since
transposes in Spot are simple book keeping operations, they are
essentially free; no underlying matrices are actually transposed.
The transpose of a complex operator, rather than its conjugate, can be
formed using the A = C.'; % Transpose of complex operator Note that in this case a single multiplication of a complex vector by
requires two multiplications by due to the implementation of
the ## Addition and subtractionThe next elementary operation is the addition and subtraction of operators. In its simplest for we add two or more operators; x = (B + C + D) * y; A = B + C + D; x = A * y; % Equivalent to first statement When Spot encounters the sum of an operator with a matrix or some
class object that implements the multiplication and size operators, it
will first wrap that entity to a Spot operator using the C = A + M; % Addition of operator and matrix C = M + A; % Addition of matrix and operator Addition of scalars to an operator is interpreted as an elementwise
addition, just like it is done for matrices. In order to make this
work we first create a new operator of appropriate size, consisting of
only ones (see A = B + 3; A = B + 3*opOnes(size(B)); Subtraction is implemented by scalar multiplication combined with addition. A = B - C; A = B + (-C); Unfortunately, Spot can only provide a limited amount of simplification of operators. This means that in the following case, no simplification is done; A = B + 3 - 2; % Not simplified, A = B + 3 - 2; A = B + C - C; % Not simplified
## Operator informationBefore we proceed with more advanced operator manipulations, let us briefly look at some ways of querying operator information. This kind of information is useful when developing algorithms and also enable us to explain how certain operator manipulations work. The most elementary property of each operator is its size, and it can
be queried using the [m,n] = size(A); % Gives m = 3, n = 6; m = size(A,1); % Gives m = 3; n = size(A,2); % Gives n = 6; o = size(A,3); % Gives o = 1; Note that the first two dimensions give the number of rows and columns
respectively. All higher dimensions have size one by definition. To
check if an operator is empty we can use if isempty(A) error('Operator A cannot be empty'); end An operator is considered empty if either its number of rows or its number of columns is zero; a operator, albeit of limited use, is perfectly valid and can be applied to vectors of corresponding size. When working interactively from the Matlab command line it is often
more convenient to display the operator by typing its name, or using
the >> A A = Spot operator: Matrix(3,6) rows: 3 complex: no cols: 6 type: Matrix >> disp(A) Spot operator: Matrix(3,6) >> whos A Name Size Bytes Class Attributes A 3x6 305 opMatrix The first two commands also provide information about the construction of the operator: B = Spot operator: Matrix(3,6)' * Matrix(3,6) rows: 6 complex: no cols: 6 type: FoG At times, for example during debugging or when using codes not
compatible with Spot, it is desirable to have access to the matrix
form underlying an operator. This can be done identically using the
>> A = opDCT(4); >> double(A) ans = 0.5000 0.5000 0.5000 0.5000 0.6533 0.2706 -0.2706 -0.6533 0.5000 -0.5000 -0.5000 0.5000 0.2706 -0.6533 0.6533 -0.2706 This explicit form is obtained though multiplication by the identity
matrix. As such, it can be quite an expensive operator. When the
number of rows is much smaller than the number of columns it may be
faster to use Finally, the command ## Dictionaries and arraysThe original motivation for developing Spot was the concatenation of operators to form a dictionary. This can now be achieved simply by writing A = [B,C,M]; % or A = [B C M]; Note that explicit matrices and classes can be mixed with Spot operators, provided they are all compatible in size. Like in operator addition Spot automatically converts these entities to Spot operators. Vertical concatenation is done likewise by typing A = [B;C;M]; % or A = [B C M]; With the specification of these two operations, Matlab automatically converts arrays of operators into a vertical concatenation of dictionaries: A = [B C; C' M]; % or A = [B C C' M]; % both represent A = [[B,C];[C',M]];
## Block diagonal operatorsAnalogous to matrices of operators it is possible to construct block
diagonal operators of operators, using the D = blkdiag(A,B,C,M); There is no restriction on the dimension of each operator. This means that the resulting operator need not necessarily be square:
> D = blkdiag(opDCT(2),ones(2,1)); double(D) ans = 0.7071 0.7071 0 0.7071 -0.7071 0 0 0 1.0000 0 0 1.0000 For more powerful constructions, including horizontal or vertical
overlap of operators the underlying command In the Operator information section we illustrated the ## Advanced interfaceThe block-diagonal operator can be invoked in two different ways: op = opBlockDiag([weight],op1,...,opn,[overlap]); op = opBlockDiag([weight],op,[overlap]); In the first mode, a block-diagonal operator is formed using operators
, with weights set to
The second mode, with only a single operator specified, can be used to
replicate the given operator. When the vectorized weight parameter,
## ExampleThe block-diagonal operator can be very useful in constructing windowed Fourier transformations.
## Kronecker productsThe D = kron(A,B,C) Needless to say, Kronecker products can increase quickly in computational complexity. Given a set of operators of size , , let , and . Then, application of the Kronecker product will require multiplications with each operator , and likewise products in transpose mode. ## ExampleKronecker products can be useful when operating on high dimensional data represented in vectorized form. For example, when represents a data volume of dimensions . To create an operator that applies a one-dimensional Fourier transformation on vectors along the second dimension, we can simply write op = kron(opEye(n),opDFT(m),opEye(l)) Likewise, the block diagonal matrix of three operators , given in
the previous section, could have been created using the Kronecker
product:
## Subset assignment and referenceFor certain applications we may be interested only in applying part of an operator. This can be done by creating new operators that are restrictions of existing operators. Restriction can be done in terms of rows, columns, a combination of rows and columns, and in terms of individual elements of the underlying matrix. Like in matrices, this is done by indexing using brackets. A = B(3:5,:); % Extract rows 3-5 A = C(:,4:6); % Extract columns 4-6 A = D(3:5,4:6); % Extract rows 3-5 of columns 4-6 The single colon ‘ A = B(3:5,1:end); % Extract rows 3-5 There is no restriction on the ordering of the selected columns and rows, and it is possible to repeat entries. A = B([1,2,2,5,3],end:-1:1); % Extract rows in reverse Because row and column indexing is implemented by pre- and post multiplication by selection matrices, the underlying operator is only applied once for each use of the restricted operator. In addition to ranges and index lists, it is also possible to use
logicals. In this mode only those entries corresponding to a
logic = randn(10,1) > 0; index = find(logic == true); A = B(logic,:); A = B(index,:); % Equivalent to previous line The orientation of the logical or index vector is irrelevant and both
row and column vectors can be used. Dimensions beyond the second one
can be added provided they are either a vector of ones, the empty set
A = B(1:3,[]); % Results in an empty 3-by-0 operator When indexing with only a single argument the numbers are interpreted
as indices in the vectorized matrix. For a given matrix,
element is accessed by absolute index . In
contrast to indexing by row and column, the result of such operations
will be a vector, rather than an operator. In fact, all required
rows or columns (depending on which one requires fewest applications
of the operator) are first computed and all relevant entries
extracted. When the only parameter is a logical it is first vectorized
and padded with x = B(:); % Vectorized entries of B y = double(B); z = y(:); % Same as x Extracting large numbers of columns or rows from an operator is computationally expensive since each column or row requires the multiplication of a vector of the identity matrix by the operator. It is thus advised to avoid using these operations whenever possible. ## ExampleAs an illustration of the use of the subset selection we create an
operator that consists of randomly selected rows of a the discrete
Fourier transform. The first step is to construct a Fourier transform
of the desired dimension using the F = opDFT(128); R = F(rand(128,1)<=0.8,:); This generates a restricted Fourier operator with approximately 80% of all rows selected. ## AssignmentExisting operators can be modified by assigning new values to a subset of the entries of the underlying matrix representation. For example we may want to zero out part of some operator : A(2,:) = 0; % Zero out row two A(:,[3,5]) = 0; % Zero out columns three and five These operations are not restricted only to zero values; other constants can be used as well. With the same ease we can also assign matrices and other operators to rectangular subsets of existing operators. In the following two examples we apply this to two matrix-based operators and , and show the resulting operator format: >> A(3:5,[2,4,6]) = randn(3,3) A = Spot operator: Subsasgn(Matrix(8,8), Matrix(3,3)) rows: 8 complex: no cols: 8 type: SubsAsgn >> B(6:-1:3,4:7) = opDCT(4) B = Spot operator: Subsasgn(Matrix(8,8), DCT(4,4)) rows: 8 complex: no cols: 8 type: SubsAsgn Note that in the second example we also flip the rows of the DCT operator upside down by choosing the row indices from 6 to 3, instead of the other way around. As long as no destination row or column is repeated (which is not permitted), and the number of selected rows and columns matches the size of the operator, any ordering can be used. The row and column indices may exceed the size of the existing
operator. Such assignments will enlarge the original operator by
logically embedding it into an operator of all zeros prior to
assigning the desired entries with the new operator, matrix or
constant (all are automatically converted to operator form). The cost
of multiplying with the resulting operator depends on whether the
embedded operator overlaps with the original operator or not. Consider
the following example (the >> A = opOnes(2,4); B = 2*opOnes(3,3); >> A(3:5,3:5) = B; >> double(A) ans = 1 1 1 1 0 1 1 1 1 0 0 0 2 2 2 0 0 2 2 2 0 0 2 2 2 In this case the operators do not overlap, and multiplication with the
new operator requires a single application of and the original
. Although the command When there is an overlap between the original and embedded operators the cost of applying the new operator requires two applications of the original and one application of , and likewise for multiplication by . >> A = opOnes(2,4); B = 2*opOnes(3,3); >> A(2:4,3:5) = B; >> double(A) ans = 1 1 1 1 0 1 1 2 2 2 0 0 2 2 2 0 0 2 2 2 In fact, multiplication by the new takes place in three stages;
the first stage performs multiplication by the original , the
second stage subtracts the influence of the overlap with , and the
third stage adds the result of multiplication by . Note that in
determining overlap Spot has to be conservative and assume that the
operators are dense. In the overlapping example above, this means that
assigning a new operator to A special situation arises when assigning the empty set >> A(:,[1,3]) = []; >> double(A) ans = 1 1 0 1 2 2 0 2 2 0 2 2 >> A([2,3],:) = []; >> double(A); ans = 1 1 0 0 2 2 While the above example uses the colon ( Assignment using absolute index numbers is not supported in the
current version of Spot. For the above operator ,
the command
## Elementwise operationsThere are a number of operations that, just like multiplication by
scalars, affect all individual entries of matrix underlying the
operator. For complex operators Spot provides the commands
F = opMatrix(2 + 3i); a = real(F); % a = opMatrix(2) b = imag(F); % b = opMatrix(3) c = conj(F); % c = opMatrix(2 - 3i) Application of
## Power, inverse, and backslash operationsThe operations discussed in this section are provided for convenience,
but they should be used with care to avoid prohibitively slow or
inaccurate operators. With that said, the power operator B = A^3 % A*A*A B = mpower(A,3) % A*A*A B = mpower(A,0) % A^0 = opOnes(size(A)) B = mpower(A,-3) % A^-3 = inv(mpower(A,3)) The B = opInverse(A); x = B * b; is computed by solving with additional regularization to obtain the least two-norm solution
in case the system is underdetermined. Because the least-squares
problem is not restricted to square matrices, The Finally, the backslash operator x = A \ b; When is underdetermined the resulting may not be the same as
## Solving systems of linear equationsMatlab provides a number of routines for solving systems of linear equations that either take a matrix or a function handle. Wrappers to these function that take Spot operators are implemented for the following functions:
## Application to non-linear operatorsIn some situations it may be desired to create nonlinear operators. For example, when the output of a certain series of operations on complex input is known to be real it may be desired to have an operator that discards the imaginary part of its input. While Spot permits the construction of such operators, utmost care has to be taken when using them in combination with meta operators. Because there are no sanity checks in the implementation of the meta operators, application to nonlinear operators is not guaranteed to be meaningful; this is something the user should decide. |