Download All

Code covered by the BSD License  

Highlights from
On-The-Fly Definition of Custom Matrix Objects

from On-The-Fly Definition of Custom Matrix Objects by Matt J
Class of matrix-like objects with on-the-fly definable methods.

UserGuide

Contents

function UserGuide

Overview of the MatrixObj and DataObj Classes.

Objects of the MatrixObj class (and its subclass DataObj) are capable of behaving as matrices, but with math operators (+,-,*,\,.*,<,>,etc...) that can be defined/redefined in any Mfile context, or even at the command line. This removes the restriction of writing a dedicated classdef file or class directory for every new matrix-type object that a user might wish to create.

The class works by storing function handles to the various matrix operator functions like plus, minus, mtimes, mldivide, etc... in a property of MatrixObj called Ops, which is a structure variable. Hence, one can set the matrix operators as desired simply by setting the fields of Ops to an appropriate function handle.

MatrixObj objects are particularly useful when an object needs to be endowed with just a few matrix-like capabilities that are very quickly expressed using anonymous functions or a few short nested functions. This is illustrated in some examples below that deal with creating an efficient version of a DFT matrix.

Another advantage of MatrixObj objects is that it is not necessary to issue a "clear classes" command when their Ops methods need to be edited or redefined. This is because the definitions of the methods are not held in the code defining the MatrixObj class. Rather, they are held in external code, to which the class is linked only through function handles.

The DataObj subclass is a specialized version of MatrixObj well-suited for mimicking/modifying the behavior of existing MATLAB numeric data types. Its Ops property contains default methods appropriate to existing data types, but which can be selectively overwritten. The section The DataObj Subclass will introduce DataObj in more detail.

MatrixObj: Preliminary Examples

We start with some basic examples using the MatrixObj class.

EXAMPLE 1: Implementing fft() in operator form.

As is well-known, the operation fft(x) can be represented as a matrix-vector multiplication. The relevant matrix can be computed as follows

 d=2500;
 Q=fft(eye(d)); %DFT matrix - 2500x2500

The operation fft(x) is equivalent to Q*x, but this is a slow way to perform the operation,

 x=rand(d);

 tic; y0=Q*x; toc

Elapsed time is 3.367615 seconds.

However, using the MatrixObj class, we can quickly create an object Qobj which can transform x using the same matrix multiplication syntax, Qobj*x, but using fft() under the hood, with all of its advantages in speed,

 Qobj=MatrixObj;
 Qobj.Ops.mtimes=@(obj,z) fft(z);  %set the mtimes method in 1 line!!


 tic; y1=Qobj*x; toc
 tic; y2=fft(x); toc

Elapsed time is 0.097657 seconds.
Elapsed time is 0.095093 seconds.

 isequal(y1,y2),   %Make sure we got the correct answer

ans =

     1

And of of course, the memory footprint of Qobj is far less than for the full matrix Q

  whos Q Qobj

  Name         Size                  Bytes  Class        Attributes

  Q         2500x2500            100000000  double       complex   
  Qobj         1x1                    4412  MatrixObj

EXAMPLE 2: Continuing with Example 1, suppose I now decide that I still want Qobj to represent an fft() operation, but that it be normalized to satisfy Parseval's theorem. A simple on-the-fly redefinition of mtimes() can accomplish this:

 Qobj.Ops.mtimes=@(obj,z) (1/sqrt(numel(z)))*fft(z);

 x=rand(d,1);

 ParsevalSatisfiedIfTheseAreEqual=[norm(x), norm(Qobj*x)],

ParsevalSatisfiedIfTheseAreEqual =

   28.7926   28.7926

EXAMPLE 3: Continuing with Example 2, let us now look at how to give Qobj a ctranspose method so that Qobj' is defined. Because Qobj satisfies Parseval's theorem, Qobj' is its inverse. There is a fair amount of flexibility in implementing (conjugate) tranposes for MatrixObj objects. This will be discussed more fully in the section Implementing Transposition. In this case, a one-line definition can be made using the Trans property, 

  Qobj.Trans.mtimes=@(obj,z) sqrt(numel(z))*ifft(z) ;

The code below verifies that the ctranpose operation has various anticipated properties,

  ParsevalSatisfiedIfTheseAreEqual=[norm(x), norm(Qobj'*x)],

  AdjointSelfCancels=isequal(Qobj*x, (Qobj')'*x),

  InversionErrorLeft=norm(x- Qobj'*(Qobj*x)),

  InversionErrorRight=norm(x- Qobj*(Qobj'*x)),

ParsevalSatisfiedIfTheseAreEqual =

   28.7926   28.7926


AdjointSelfCancels =

     1


InversionErrorLeft =

  8.9156e-015


InversionErrorRight =

  9.0291e-015

MatrixObj: Details of Usage

The following constructs a default MatrixObj object M,

    M=MatrixObj;

The two main properties of MatrixObj are "Params" and "Ops". The properties can be populated with non-default values by dot-indexing assignment statements. As shown in the examples above, there is also a "Trans" property, which will be discussed later.

Another construction syntax is

   M=MatrixObj(Params,Ops)

which can be used to set these properties directly.

The Params property can be any MATLAB variable - numeric, struct, or otherwise - and is meant for holding object parameters of any kind. Its default/initial value is [].

The Ops property is a structure whose field names include all of the math operator functions, plus a few more such as as sum(), inv(), conj(), size(), and display:

   and             horzcat       mpower        power          sum
   colon           inv           mrdivide      rdivide        times
   ctranspose      ldivide       mtimes        size           transpose
   end             le            ne            subsasgn       uminus
   eq              lt            not           subsindex      uplus
   ge              minus         or            subsref        vertcat
   gt              mldivide      plus          display        conj

By assigning appropriate function handles to these fields, any/all of the corresponding class methods can be implemented on-the-fly. The fields of Ops also default to [].

Although M.Ops is in most ways manipulable as a structure using dot-indexing syntax, any attempt to assign a field to Ops which is not in methods(MatrixObj) will result in an error. However, advanced users, if they wish, can add more methods to the MatrixObj classdef, and these methods will then automatically join the list of manipulable fields of Ops.

Implementing Transposition

In many applications, meaningful transposition of a MatrixObj requires extensive redefinition of the object and its operations. To facilitate this, the ctranspose operation M' and the transpose M.' are implemented with some special handling and features. In particular, these operations trigger the following steps:

(1) First, M.Ops.ctranspose is executed (or Ops.tranpose, as appropriate) on M.

(2) If the result of (1) is also a MatrixObj, its Trans and Ops properties are swapped. The Trans property, when set non-empty, is a struct of the same format as Ops, i.e., it has the same field names. It is meant to specify an alternative set of math operations for the (conjugate) tranpose of M.

(3) As explained in the section More Special-Case MatrixObj Methods, one has the option of using a numeric variable [p,q] to specify Ops.size of a MatrixObj instead of a function handle. If M.Ops.size holds a numeric variable [p,q], then p and q will be swapped so that M' and M.' will have Ops.size=[q,p].

If Ops.ctranspose=[] as is the default, step (1) is skipped, leaving M unaltered. Similarly, if Trans=[], the swapping in step (2) is omitted. Therefore, one has the flexibility of using either, both, or neither of these two class properties to contribute to the implementation of a ctranspose method. When both Ops.ctranspose=[] and Trans=[], neither will be used, and no change in the object will result from the M' operation. Thus, all MatrixObj objects are "symmetric" and "conjugate symmetric" by default. The dimension swapping in step (3) will only ever be executed if M.Ops.size is numeric.

In Example 3, only the Trans property was used to implement the ctranspose. The implementation consists entirely of step (2). To see its effect, we continue from Example 3 and create

 Qtobj=Qobj';

We can now verify directly that Qtobj.Ops.mtimes and Qtobj.Trans.mtimes are the same as for Qobj, but interchanged.

 OpsMtimes = Qtobj.Ops.mtimes,
 TransMtimes = Qtobj.Trans.mtimes,

OpsMtimes = 

    @(obj,z)sqrt(numel(z))*ifft(z)


TransMtimes = 

    @(obj,z)(1/sqrt(numel(z)))*fft(z)

The Trans property sometimes gives a conveniently quick way of implementing a (c)transpose. As Example 3 showed, a conjugate tranpose transformation was implemented this way in just one line of code. However, it doesn't always offer the desired flexibility. The next Example shows a situation where the Ops.ctranpose property provides a better alternative.

EXAMPLE 4: In this example, we create a MatrixObj used to imitate a diagonal matrix, but which internally stores the diagonal elements only. This example implements the conjugate transpose using Ops.ctranspose exclusively.

  M=MatrixObj;

   M.Params=[1:5]*i; %Use the Params property to hold the diagonal data
                     %This will be the same as diag(1:5)*i;

   M.Ops.mtimes=@(obj,x) (obj.Params(:)).*x; %Multiplication method
   M.Ops.ctranspose=@(obj) set(obj,'Params',conj(obj.Params));

Now let's try a few operations

  MultOp = M*(1:5)'
  AdjointMultOp = M'*(1:5)'

MultOp =

        0 + 1.0000i
        0 + 4.0000i
        0 + 9.0000i
        0 +16.0000i
        0 +25.0000i


AdjointMultOp =

        0 - 1.0000i
        0 - 4.0000i
        0 - 9.0000i
        0 -16.0000i
        0 -25.0000i

EXAMPLE 5: This additional example is a continuation of the FFT operator application. Here, we use Ops.ctranpose and Trans in concert to implement the conjugate tranpose. The implementation will be the same as in Example 3, except that now, we will use the Params property to label whether the operator is a forward or inverse FFT. The forward FFT will be labeled with +1 and the inverse will be labeled -1. Although an abstract example for now, it will serve a more meaningful purpose in later examples.

 Qobj=MatrixObj;
 Qobj.Params=+1; %the label

 Qobj.Ops.mtimes=@(obj,z) (1/sqrt(numel(z)))*fft(z); % still Parseval normalized
 Qobj.Trans.mtimes=@(obj,z) sqrt(numel(z))*ifft(z) ;

 minusParams=@(obj) set(obj,'Params',-obj.Params);
 Qobj.Ops.ctranspose=minusParams;
 Qobj.Trans.ctranspose=minusParams;

Thus, the change in the Params label is handled by Ops.ctranpose, whereas the implementation of the actual transform (the swap from fft to ifft) is implemented through the Trans property. The next few test operations illustrate this.

 Qtobj=Qobj'; %perform a ctranpose operation

 fftMtimes = Qobj.Ops.mtimes, %verify change in the mtimes implementation
 ifftMtimes = Qtobj.Ops.mtimes,

 fftLabel = Qobj.Params, %verify change in label
 ifftLabel =Qtobj.Params,

fftMtimes = 

    @(obj,z)(1/sqrt(numel(z)))*fft(z)


ifftMtimes = 

    @(obj,z)sqrt(numel(z))*ifft(z)


fftLabel =

     1


ifftLabel =

    -1

Over-riding subsref and subsasgn for MatrixObj

By default, the only indexing operations enabled for the MatrixObj class is dot-indexing, and is provided as a way of accessing properties. If desired however, Ops.subsref and Ops.subsasgn can be used to redefine indexing operations, for example if one desires to implement () and {} indexing as well. However, there are some subtle difficulties with this that require some care. This is much less a problem for the DataObj subclass, see Over-riding subsref and subsasgn for DataObj

Firstly, if M.Ops.subsref/subsasgn are assigned handles, the default dot-indexing capability will no longer be available for property access commands like "M.Ops" and "M.Params". If it is desired to continue to have dot-indexing access to properties, M.Ops.subsref and/or M.Ops.subsasgn must reimplement dot-indexing appropriately. Otherwise, the MatrixObj class has set() and get() methods which can always be used to access the properties of M.

Secondly, it is essential that the functions pointed to by the handles Ops.subsref/subsasgn do not contain any indexing expressions of the MatrixObj object. Any such expressions will cause an infinite recursive call to Ops.subsref/subsasgn. Thus, MatrixObj properties must always be accessed using set() and get() from within the functions used to redefine subsref/subsasgn.

More Special-Case MatrixObj Methods

By default, most fields of Ops start with Ops.methodname=[] meaning that the corresponding method is unimplemented. If called, unimplemented methods will normally return an error. In the previous sections, we saw however that Ops.transpose(), Ops.ctranpose(), Ops.subsref(), and Ops.subasgn() are special cases with specific default behavior. There are a few more special case methods to be aware of:

  • Ops.display - When empty, the MATLAB's default display of the object will be used.
  • Ops.conj - When empty, conj(M) will leave MatrixObj M unchanged. Hence, every MatrixObj is "real" by default.
  • Ops.size - When empty, size(M) will return [1,1] by default. Hence, every default MatrixObj is a "scalar". This is also the one method that can be specified as a numerical variable instead of a function handle. If M.Ops.size=[p,q], then size(M) returns [p,q]. In this case also, transpose operations will give the Ops.size field special handling, as described earlier. Finally, it is never necessary to implement multiple input arguments for Ops.size. As long as size(M) is implemented to return a [p,q], whether by function handle or by storing [p,q] explicitly, then size(M,1) returns p and size(M,2) returns q.
  • Ops.end - When empty, expressions like M(:,end) will attempt to derive the value of "end" from size(M). Of course, this assumes that Ops.subsref or similar has been used to supply a definition for ()-indexing.

EXAMPLE 6: In this example, we add size() and subsref() methods to the FFT operator defined in Example 5 to illustrate some of the above special cases. To implement a subsref method, we will use a nested function

Qobj.Ops.size=[5 5];  %We'll use a smaller dimension than before
Qobj.Ops.subsref=@NestedSubsref;
Qobj.Trans.subsref=@NestedSubsref;


    function out=NestedSubsref(obj,S)

        switch S.type

            case '.'  %re-implement dot-indexing notation

               out=get(obj,S.subs); %I can't execute this as obj.(S.subs)
                                    %or it will cause infinite recursion.
                                    %I have to use the class get() method!!!

            case '()'

              ExpSign=get(obj, 'Params'); %Either +1 or -1 depending on FFT or IFFT

              N=size(obj,1);

              m=(1:N).';      %parse subscripts obj(m,n) - we will not implement linear indexing
               m=m(S.subs{1});
              n=1:N;
               n=n(S.subs{2});

              out=1/sqrt(N)*exp(-ExpSign*j*2*pi/N*(m-1)*(n-1));

        end

    end

Now, let's try some operations

Qtobj=Qobj';

Dimensions = size(Qobj),
Dimension1 = size(Qobj,1),
FinalElement = Qobj(end,end),

ColumnsShouldBeEqual = [Qobj(:,3), Qobj*[0;0;1;0;0]],

ColumnsShouldBeConjugates = [Qobj(:,3), Qtobj(:,3)],

ExpSign = Qobj.Params, %Check to see that dot-indexing still available

Dimensions =

     5     5


Dimension1 =

     5


FinalElement =

   0.1382 - 0.4253i


ColumnsShouldBeEqual =

   0.4472             0.4472          
  -0.3618 - 0.2629i  -0.3618 - 0.2629i
   0.1382 + 0.4253i   0.1382 + 0.4253i
   0.1382 - 0.4253i   0.1382 - 0.4253i
  -0.3618 + 0.2629i  -0.3618 + 0.2629i


ColumnsShouldBeConjugates =

   0.4472             0.4472          
  -0.3618 - 0.2629i  -0.3618 + 0.2629i
   0.1382 + 0.4253i   0.1382 - 0.4253i
   0.1382 - 0.4253i   0.1382 + 0.4253i
  -0.3618 + 0.2629i  -0.3618 - 0.2629i


ExpSign =

     1

The DataObj Subclass

Sometimes it is desirable to customize or modify the behavior of other existing data types or classes. The DataObj subclass is tailored to situations like these. Instead of a "Params" property, this subclass contains a "Data" property for holding an object data variable. Typically, this would be an existing MATLAB numeric variable type, although this is not strictly necessary. When a function Ops.MethodName is called, the function will receive this Data variable as input, in place of the actual DataObj object. This is in contrast to MatrixObj Ops functions, which always receive the entire object as input. The result of a DataObj operation is also always a new DataObj object, one whose Data property holds the results of the computation. The function Ops.size() is an exception to this rule. The size() method will return normal numeric output. These ideas are demonstrated in the following example.

EXAMPLE 7: In this example, DataObj is used to create a specialized array type which invokes bsxfun() for a variety of operations. This can be a useful way of circumventing bsxfun's lengthy functional syntax.

 P=DataObj;
  P.Data=[1,2;3,4].',

        P.Ops.minus=@(A,B) bsxfun(@minus,A,B);
        P.Ops.plus= @(A,B) bsxfun(@plus,A,B);
        P.Ops.times=@(A,B) bsxfun(@times,A,B);
        P.Ops.rdivide= @(A,B) bsxfun(@rdivide,A,B);
        P.Ops.ldivide= @(A,B) bsxfun(@ldivide,A,B);
        P.Ops.power=@(A,B) bsxfun(@power,A,B);
        P.Ops.eq=@(A,B) bsxfun(@eq,A,B);
        P.Ops.ne= @(A,B) bsxfun(@ne,A,B);
        P.Ops.lt=@(A,B) bsxfun(@lt,A,B);
        P.Ops.le= @(A,B) bsxfun(@le,A,B);
        P.Ops.gt= @(A,B) bsxfun(@gt,A,B);
        P.Ops.ge=@(A,B) bsxfun(@ge,A,B);
        P.Ops.and= @(A,B) bsxfun(@and,A,B);
        P.Ops.or=@(A,B) bsxfun(@or,A,B);


  Q=P-[1,2],
  R=P+[3;7],

P =

     1     3
     2     4
 

Q =

     0     1
     1     2
 

R =

     4     6
     9    11

Notice that the above P.Ops functions are written in terms of normal MATLAB numerical input variables. There is no need for these functions to extract the Data to be operated upon from the DataObj. This extraction is done by the DataObj class methods automatically. Contrast this with Example 4, in which each M.Ops function was written in terms of MatrixObj input arguments. In that example, the supplied M.Ops functions therefore had to dig the Params values out of the object first before operating on them.

Notice also that the Q and R objects resulting from these operations display on the screen like normal MATLAB variables. However, they are in fact DataObj objects as this shows,

  whos P Q R

  Name      Size            Bytes  Class      Attributes

  P         2x2              4908  DataObj              
  Q         2x2              4908  DataObj              
  R         2x2              4908  DataObj

and they have the same bsxfun-driven methods.

  Q.Ops.minus,
  R.Ops.plus,

ans = 

    @(A,B)bsxfun(@minus,A,B)


ans = 

    @(A,B)bsxfun(@plus,A,B)

In general, operations on DataObj method, except for size(), will result in a copy of the object that dispatched the method, but with altered Data.

Another difference between DataObj and MatrixObj is that in MatrixObj, all members Ops.MethodName default to [], so that most methods will fail to execute when an implementation is not explicitly provided. Conversely, DataObj objects get constructed with natural defaults. For instance, even though we did not specify an mtimes implementation for P, the usual mtimes is available by default and operates on P.Data in the natural way.

   Z=P*diag(1:2),

Z =

     1     6
     2     8

This is also why the objects P, Q, and R display to the screen like normal MATLAB variables. It is the natural display() default at work. What is in fact being displayed are the Data properties of these objects,

P.Data, Q.Data, R.Data,

ans =

     1     3
     2     4


ans =

     0     1
     1     2


ans =

     4     6
     9    11

EXAMPLE 8: In this example, motivated by a newsgroup discussion, we create a DataObj which is in most ways like a normal numeric MATLAB variable, except that comparisons A==B return true if they are satisfied within a tolerance of TOL=.001.

 P=DataObj;
  P.Data=2;
  TOL=.001;

  P.Ops.eq=@(A,B) abs(A-B)<=TOL;

  P==2,        %true
  P==2.0001,   %true
  P-1==1.0001, %true
  P==2.01,     %false

ans =

     1
 

ans =

     1
 

ans =

     1
 

ans =

     0

Over-riding subsref and subsasgn for DataObj

Unlike MatrixObj, dot-indexing is always available in the DataObj subclass, and is not redefinable via Ops.subsref and Ops.subsasgn. It can be used for class property access only. This restriction prevents DataObj from being used to customize struct data types. DataObj is intended primarily for customizing numeric types.

Also in contrast to MatrixObj, ()-indexing and {}-indexing have a default implementation for any DataObj object, D. These kinds of indexing expressions will call the subsref/subsasgn method of D.Data, if such a method exists.

Finally, when index expressions D(...) or D{...} are redefined through function handles Ops.subsref or Ops.subsasgn, the functions will be passed D.Data, just like all the other methods. There is no need to use get() and set() to extract the D.Data in these implementations.

More Special-Case DataObj Methods

For MatrixObj objects, the following methods had special behavior and defaults. For a DataObj object, D, these defaults get revised as follows,

  • D.Ops.display - When empty, the class will call display(D.Data).
  • D.Ops.conj - Will attempt to implement as conj(D.Data).
  • D.Ops.size - Will attempt to implement as size(D.Data). The output is always numeric, unlike other DataObj methods.
  • Ops.end - Same default behavior as for MatrixObj. If Ops.end is a non-empty handle, D.Data will be passed to D.Ops.end.
  • Ops.transpose, Ops.ctranspose - Will try to invoke tranpose(D.Data) or ctranspose(D.Data) as appropriate. If Trans property is non-empty, swapping will occur as with MatrixObj.
end


Published with MATLAB® 7.9

Contact us at files@mathworks.com