ARRAYS

In addition to the basic ${\tt <type\gt space}$ and ${\tt <type\gt op}$ classes, there are also classes that define arrays of each. A ``${\tt <type\gt spacearray}$'' is a two dimensional array of ${\tt <type\gt spaces}$, which are not necessarily conformable with each other. An operator array is also a two dimensional array which contains operators. Applying a ${\tt <type\gt oparray}$ to a ${\tt <type\gt spacearray}$:

$${\bf u} = {\bf A v}$$

works like normal matrix multiplication. All three of ${\bf u}$, ${\bf v}$, and ${\bf A}$ are two dimensional arrays (of operators or spaces.)

$$\pmatrix{ {\bf u}_{11} & \cdots & {\bf u}_{1n} \cr \vdots &... ... \vdots & \ddots & \cr {\bf v}_{k1} & \cdots & {\bf v}_{kn} }$$

u is a space array that has the number of rows of A, and the number of columns of v. It is not possible to apply an operator array to an array of unsuitable size without generating a run time error. The rows of A are applied to the columns of v and the result is summed. As each element of the operator array is applied to the appropriate element of the space array, run time compatibility checks will take place. e.g. ${\bf A}_{11}$ and ${\bf v}_{11}$ must be compatible. Also the results to be summed must be compatible, the space that results from ${\bf A}_{11} {\bf v}_{11}$must be compatible with the space that results from ${\bf A}_{12} {\bf v}_{21}$.

In C++, this looks like:

<type>spacearray u = B.Forward(v);

This is almost identical to the syntax for applying a single operator to a single space, except here the variables u and v are of class ${\tt <type\gt spacearray}$ and the operator B is of class ${\tt <type\gt oparray}$.

A simple example of using arrays is to generate a damped operator from the regular operator. If we wish to calculate the damped solution to the problem

$$\min_x \left\vert{\bf A x} - {\bf y} \right\vert^2$$

We can augment the problem with a diagonal damping matrix and solve:  

$$\min_x \left\vert \pmatrix{{\bf A} \cr \epsilon {\bf I} } {\bf x} - \pmatrix{ {\bf y }\cr {\bf 0} } \right\vert^2$$ (1)

The new operator, ${\bf A_2}$ can be written as

$${\bf A_2} = \pmatrix{ {\bf A} \cr \epsilon {\bf I} }$$

When applied to a space, ${\bf x}$ it gives a space-array,

$${\bf A_2} {\bf x} = \pmatrix{ {\bf Ax} \cr \epsilon {\bf I x } }$$

When the conjugate is applied to a space-array it gives a space,

$${\bf A_2^{\dagger} } \pmatrix{ {\bf y_1} \cr {\bf y_2} } = ... ...{\bf y_2} } = {\bf A^\dagger y_1 }+ \epsilon^* {\bf I}{\bf y_2}$$

All of this is handled automatically by the ${\tt <type\gt oparray}$ class.

The following routine gives the damped least squares solution to any linear problem that can be posed in terms of opa, an instance of the class ${\tt <type\gt op}$,and a RHS, y, of the class ${\tt <type\gt space}$.

// This routine gives a damped solution to the least squares problem 
// opa.Forward( x ) = y;
// by augmenting the operator with a diagonal damping matrix 

<type>space damped_solver( <type>space& x0,   // initial solution
                           <type>op opa,      // linear operator
                           <type>space& y,    // rhs
                           float epsilon,     // the damping factor 
                           int maxiter )      // maximum no. of iterations
{
  <type>oparray ops(2);            // declare a 2x1 operator array
  <type>Scale   scale( epsilon );  // scale is a diagonal scaling operator
  ops(1) = &opa; ops(2) = &scale; // assign the elements of the operator array

  <type>spacearray rhs(2)          // declare a 2x1 rhs array of spaces.
  <type>space zero = x.empty()     // zero is the same shape as x but all zeros.
  rhs(1) = y; rhs(2)=zero         // assign the elements of the rhs

  // now solve the augmented problem, and return the result
  return <type>space soln = <type>arraysolver( x0, ops, rhs, maxiter );

}

The new operator ops ( of class ${\tt <type\gt oparray}$) and right hand side, rhs ( of class ${\tt <type\gt spacearray}$ ) are used to solve the damped least squares problem. This is done by the routine ${\tt <type\gt arraysolver}$, which is one of a family of solver routines that are available.