MatrixLib  

Site Index:

MatrixLib

Contents

1. Introduction

Back to MatrixLib Table of Contents   OptiVec home

1.1 Matrix Data Types

The ordering of elements is the same as in the two-dimensional arrays provided by the respective target compilers. This means that the matrices are stored row-wise in MatrixLib versions for C and C++ compilers, but column-wise in the Pascal/Delphi versions.

While we recommend to exclusively use these dynamically allocated matrix types, static matrices defined, e.g., as
float MX[4][6]; (for C/C++), or
MX: array[0..3][0..5] of Single; (for Pascal/Delphi)
can be used in all MatrixLib functions with the exception of the multiLinfit and multiNonlinfit routines.

Back to MatrixLib Table of Contents   OptiVec home

1.2 Prefixes of MatrixLib Functions

• MFdia_ means that the function expects a diagonal matrix (i.e., a square matrix which has non-zero elements only on the diagonal); as there is no sense in storing all the zeros, diagonal matrix are actually stored as vectors, holding only the diagonal elements.
   
• MFsym_ denotes a function which expects the input matrix to be symmetric. At present, only MFsym_eigenvalues makes use of this assumption.
   
• MFtrd_ means the function is for a tridiagonal matrix (i.e., a square matrix with non-zero elements only on the diagonal plus or minus one column). A tridiagonal matrix has to be entered in the form of a matrix with three rows, representing the three vectors actually containing non-zero data. In other words, the original matrix
   

  d0	u0	0	···
  l1	d1	u1	0	···
  0	l2	d2	u2	0	···
  			···
  			lN-2	dN-2	uN-2
  			0	lN-1	dN-1

   
  is compacted into the form
   

  u0	u1	u2	···	uN-2	*
  d0	d1	d2	···	dN-2	dN-1
  *	l1	l2	···	lN-2	lN-1

   
  The elements l₀ and u_N-1, marked with an asterisk, are undefined and never used.

Back to MatrixLib Table of Contents   OptiVec home

1.3 C/C++ Version Specifics:

1.3.1 MatObj, the object-oriented interface for matrix functions

By analogy with the alias names of the vector objects, the matrix objects get the alias names fMatObj, dMatObj, and so on, with the data-type signalled by the first one or two letters of the class name.

The matrix constructors are:
matrix(); // no memory allocated, dimensions set to 0
matrix( unsigned ht, unsigned len ); // ht by len matrix allocated
matrix( unsigned ht, unsigned len, T fill ); // as before, but initialized with value "fill"
matrix( matrix <T> init ); // creates a copy of the matrix "init"

For the matrix classes, the arithmetic operators
+    -    *    +=    -=    *=
are defined.
While the vector-vector operator * refers to element-wise multiplication, the matrix-matrix, vector-matrix, or matrix-vector operator * means true matrix-matrix or matrix-vector multiplication.

If you ever need to process a MatObj matrix in a "classic" plain-C OptiVec function, you may use the member functions
getHt() and getLen() to retrieve the matrix dimensions,
getMatrix to obtain the actual matrix pointer of the data type fMatrix etc., or
getM0 for the vector M[0] of the data type fVector etc., which is actually a pointer to the first matrix element.

The syntax of all MatObj functions is described below together with the basic MatrixLib functions for which tMatObj serves as a wrapper.

Please note the following restrictions, coming from the tight control of matrix dimensions for compatibility:

• With the "normal" MF_ functions, one can use some large matrix as working space for matrices of smaller dimensions. As long as one does not access individual elements or rows, but uses only whole matrices, this allows to use one and the same working space, say for an [n * n] and an [m * n] matrix (m < n: neglect one or more rows). This is not possible in MatObj, where the matrix dimensions are encapsulated, and no distinction between "actual" and "used" dimensions can be made.
• It is not even possible to overwrite, e.g., a non-square matrix with its transpose, even if you were willing to sacrifice the row pointers, and even if the number of matrix elements obviously stays the same.
• All MatObj Dia_... functions require square matrices (ht = len), whereas MF_Dia_... functions work, in principle, with non-square matrices, as long as the parameter "len" in the function call corresponds to the smaller dimension.
• The multifit functions are not encapsulated, as the VF_EXPERIMENT and MF_EXPERIMENT structures need pointers rather than objects anyway.

1.3.2 Direct-pointer form of matrix function calls

Back to MatrixLib Table of Contents   OptiVec home

1.4 Pascal/Delphi Version Specifics:

The situation is somewhat different for dynamic arrays in Delphi. While, in the one-dimensional case, they they can be used with all VectorLib functions by simply type-casting fArrays into fVectors, two-dimensional arrays require a different approach: Dynamic two-dimensional arrays can be created in Delphi by declaring them as "array of fArray", "array of dArray", etc. Short-cut definitions for these types are also given in the VecLib unit, as f2DArray, d2DArray, etc. As a consequence of the described way of defining 2D-Arrays in Delphi, each row of a matrix is stored in a separate vector. This means that the matrix elements do no longer occupy a single, contiguous memory space. Therefore, 2DArrays cannot be used directly with MatrixLib functions. Instead, they have to be copied into matrices by calling MF_2DArrayToMatrix etc., before MatrixLib functions can be used on them. The reverse conversion is available as MF_MatrixTo2DArray.

In the following chapters, a brief summary of all MatrixLib function is given, ordered into groups according to their functionality. Click on the links to obtain more information about the individual functions mentioned.

Back to MatrixLib Table of Contents   OptiVec home

2. Management of Dynamically Generated Matrices

C/C++ version only:
OptiVec's dynamically allocated matrices can be addressed just like two-dimensional static arrays of C/C++. If you have, e.g., an fMatrix MX; and a variable float a;, you can write a line like
a = MX[3][5];

Both C/C++ and Pascal/Delphi versions:
There are two functions addressing single elements. These functions are the only means to address elements in Pascal/Delphi, and are also needed for getting around the pointer arithmetics bug in some versions of Borland C++:

MF_Pelement	Pointer to a specific matrix element
MF_element	value of a specific matrix element

To assign a value to a specific matrix element, you should use the syntax
*(MF_Pelement( MX, ht, len, 3, 4 )) = 3.0; /* for C/C++*/
MF_Pelement( MX, ht, len, 3, 4 )^ := 3.0; (* for Pascal/Delphi *)

Back to MatrixLib Table of Contents   OptiVec home

3. Initialization of Matrices

Back to MatrixLib Table of Contents   OptiVec home

4. Data-Type Conversions

Back to MatrixLib Table of Contents   OptiVec home

5. Symmetry Operations (Transposing, Rotating, Reflecting),
Interpolation, Augmenting, Deleting, Extracting, and Filling Parts of a Matrix

Back to MatrixLib Table of Contents   OptiVec home

6. Arithmetic Operations Performed on a Single Row, Column, or the Diagonal

Back to MatrixLib Table of Contents   OptiVec home

7. Operations Performed Along All Rows or All Columns Simultaneously, or Along the Diagonal of a Square Matrix
Calculation of the Center of Gravity

Please note that multiplying all rows or all columns by one and the same vector is equivalent to a multiplication by a diagonal matrix, which is provided by MF_mulMdia and MF_TmulMdia.
For all of the above functions involving maxima (...max,  ...absmax), the corresponding minima are found by the functions named ...min,  ...absmin.
For complex numbers, maxima can be defined by various criteria. The following table summarizes the functions finding them (again, of course, the corresponding functions for the minima exist as well):
 

MCF_Rows_absmax	store the maxima of the absolute values of all rows in a (real-valued) column vector
MCF_Cols_absmax	store the maxima of the absolute values of all colums in a (real-valued) row vector
MCF_Dia_absmax	return the maximum of the absolute values of the diagonal elements as a (real) scalar
MCF_Rows_absmaxReIm	find the maximum absolute values of all real and of all imaginary parts separately along the rows of a matrix; merge the maxima into complex numbers, and store them in a column vector
MCF_Cols_absmaxReIm	find the maximum absolute values of all real and of all imaginary parts separately along the columns of a matrix; merge the maxima into complex numbers, and store them in a row vector
MCF_Dia_absmaxReIm	find the maximum absolute values of all real and of all imaginary parts separately along the diagonal of a square matrix; merge these two maxima into one complex number and return it as a scalar
MCF_Rows_cabsmax	find the complex numbers of largest magnitude along rows and store these row maxima in a (complex) column vector
MCF_Cols_cabsmax	find the complex numbers of largest magnitude along columns and store these column maxima in a (complex) row vector
MCF_Dia_cabsmax	find the complex number of largest magnitude along the diagonal of a square matrix and return it as a (complex) scalar
MCF_Rows_maxReIm	find the maxima of all real and of all imaginary parts separately along the rows of a matrix; merge the maxima into complex numbers, and store them in a column vector
MCF_Cols_maxReIm	find the maxima of all real and of all imaginary parts separately along the columns of a matrix; merge the maxima into complex numbers, and store them in a row vector
MCF_Dia_maxReIm	find the maximum of all real and of all imaginary parts separately along the diagonal of a square matrix; merge these two maxima into one complex number and return it as a scalar
MCF_Rows_sabsmax	find the complex numbers of largest sum |Re|+|Im| along rows and store these row maxima in a (complex) column vector
MCF_Cols_sabsmax	find the complex numbers of largest sum |Re|+|Im| along columns and store these column maxima in a (complex) row vector
MCF_Dia_sabsmax	find the complex number of largest sum |Re|+|Im| along the diagonal of a square matrix and return it as a (complex) scalar

To determine the center of gravity of a matrix, you have the choice between the following two functions:
 

MF_centerOfGravityInd	center of gravity, returned as interpolated element indices
MF_centerOfGravityV	center of gravity of an MZ matrix with explicitly given X and Y axes

Back to MatrixLib Table of Contents   OptiVec home

8. Operations Involving Two Rows or Two Colums

Back to MatrixLib Table of Contents   OptiVec home

9. Whole-Matrix Arithmetics: Addition, Multiplication

b)Matrix multiplication:

MF_mulV	multiply a matrix by a column vector: Y = MA * X
MF_TmulV	multiply the transpose of a matrix by a column vector: Y = MAT * X
VF_mulM	multiply a row vector by a matrix: Y = X * MA
VF_mulMT	multiply a row vector by the transpose of a matrix: Y = X * MAT
MF_mulM	multiply two matrices: MC = MA * MB
MF_mulMT	multiply one matrix by the transpose of another matrix: MC = MA * MBT
MF_TmulM	multiply the transpose of a matrix by another matrix: MC = MAT * MB
MF_TmulMT	multiply the transpose of one matrix by the transpose of another matrix: MC = MAT * MBT
MCF_mulMH	multiply one matrix by the hermitian conjugate of another matrix: MC = MA * MBT *
MCF_HmulM	multiply the hermitian conjugate of a matrix by another matrix: MC = MAT * * MB

c) Multiplications of general matrices by diagonal matrices or vice versa
The diagonal matrix is passed to the respective function as a vector.

MF_mulMdia	multiply a general matrix by a diagonal matrix: MC = MA * MBDia
MF_TmulMdia	multiply the transpose of a matrix by a diagonal matrix: MC = MAT * MBDia
MFdia_mulM	multiply a diagonal matrix by a general matrix: MC = MADia * MB
MFdia_mulMT	multiply a diagonal matrix by the transpose of a general matrix: MC = MADia * MBT

Back to MatrixLib Table of Contents   OptiVec home

10. Linear Algebra

Back to MatrixLib Table of Contents   OptiVec home

11. Eigenvalues and Eigenvectors, Matrix Square-Root

Back to MatrixLib Table of Contents   OptiVec home

12. Two-Dimensional Fourier-Transform Methods

Back to MatrixLib Table of Contents   OptiVec home

13. Data Fitting

13.1 Polynomials

13.2 General Linear Model Functions

void LinModel( fVector BasFuncs, float x, unsigned nfuncs )
{
  BasFuncs[0] = sin(x);
  BasFuncs[1] = cos(x);
}

Note that the coefficients a_i are not used in the model function itself. The argument nfuncs (which is neglected in the above example) allows to use a variable number of basis functions. It is also possible to switch fitting parameters "on" and "off", i.e., "free" or "fixed". To this end, the routine VF_linfit takes a "status" array as an argument. For all members of the parameter vector that are to be treated as free, the corresponding "status" entry must be set to 1. If status_i is set to 0, the corresponding parameter, a_i, is treated as fixed at its value upon input.
Internally, VF_linfit employs a Singular Value Decomposition algorithm to obtain a solution even for (near-)singular linear systems. Thereby, coefficients a_i whose significance is lower than a threshold, Thresh, are set to 0 instead of infinity. This threshold can be modified by calling VF_setLinfitNeglect. The current threshold can be retrieved by VF_getLinfitNeglect.

Having noted above that functions like
y = sin(a₀x)+ cos(a₁x)
require their own treatment, we arrive at the last and most general class of model functions:

13.3 Non-Linear Models

• the Levenberg-Marquardt method, and the
• Downhill-Simplex method by Nelder and Mead.

void NonlinModel( fVector Y, fVector X, ui size, fVector A)
{
  for( ui i=0; i<size; i++ )
  Y[i] = sin( A[0]*X[i] ) + cos( A[1]*X[i] );
}

VF_nonlinfit will call this model function again and again, hundreds and thousands of times, with different values in the parameter vector A. This means that a nonlinear fitting routine will spend most of its time in your model function (and its derivatives, see below), and you might wish to optimize this model function as far as possible, using vector functions. Suppose you have declared and allocated the fVector YHelp somewhere within your main program. Then you can write:

void OptNonlinModel( fVector Y, fVector X, ui size, fVector A)
{
  VFx_sin( YHelp, X, size, A[0], 0.0, 1.0 );
  VFx_cos( Y,     X, size, A[1], 0.0, 1.0 );
  VF_addV( Y, Y, YHelp, size );
}

The parameter array A you just saw being used by VF_nonlinfit to call your model function, is the one you have to pass as the first argument to VF_nonlinfit. Upon input, "A" must (!) contain your initial "best guess" of the parameters. The better your guess, the faster VF_nonlinfit will converge. If your initial guess is too far off, convergence might be very slow or even never attained. Upon output, "A" will contain the best set of parameters that could be found.

All nonlinfit functions return the figure-of-merit of the best parameter array "A" found. To perform the fit, these functions need not only the fitting function itself, but also its partial derivatives with respect to the individual coefficients. Therefore, an argument "Derivatives" is required in calling the nonlinfit functions. If you happen to know the partial derivatives analytically, Derivatives should point to a function calculating them. If you do not know them, call with Derivatives = NULL (nil in Pascal/Delphi). In the latter case, all derivatives will be calculated numerically inside the functions. In cases where you know some, but not all of the partial derivatives of Y with respect to the coefficients a_i, you could also calculate dY / da_i wherever you have an analytic formula, and call VF_nonlinfit_autoDeriv for the remaining coefficients. To demonstrate this possibility, here is a function coding the derivatives for our non-linear sine and cosine example:

void DerivModel( fVector dYdAi, fVector X, ui size, unsigned iPar, fVector A, VF_NONLINFITWORKSPACE *ws )
{
  switch( iPar )
  {
    case 0:
      for( ui i=0; i<size; i++ )
        dYdAi[i] = X[i] * cos( A[0]*X[i] );
      break;
    case 1: /* say we don't know this derivative: */
      VF_nonlinfit_autoDeriv( dYdAi, X, size, ipar, A, &ws );
  }
}

As you might have guessed, the only purpose of the input argument "ws" for DerivModel is to be passed on to VF_nonlinfit_autoDeriv, in case this becomes necessary. Think twice, however, before resorting to VF_nonlinfit_autoDeriv: The function calls needed to determine the derivatives numerically may easily slow down your fit by a factor of 3 to 10! Anyway, as described above for the model function itself, also its derivatives should be implemented using vector functions. This leads us finally to the optimized version of the derivatives function:

void OptDerivModel( fVector dYdAi, fVector X, ui size, unsigned iPar, fVector A, VF_NONLINFITWORKSPACE *ws )
{
  switch( iPar )
  {
    case 0:
      VFx_cos( dYdAi, X, size, A[0], 0.0, 1.0 );
      VF_mulV( dYdAi, dYdAi, X, size ); break;
    case 1:
      VFx_sin( dYdAi, X, size, A[1], 0.0, -1.0 );
      VF_mulV( dYdAi, dYdAi, X, size );
  }
}

Actual working examples for the polynomial, general linear, and general non-linear cases are contained in the OptiVec shareware versions you can download.

The nonlinear fitting routines are highly sophisticated and offer the user a lot of different options. These options may be set by the function V_setNonlinfitOptions. To retrieve current settings, use
V_getNonlinfitOptions.
All options are packed into a structure (struct in C/C++, record in Pascal/Delphi) named VF_NONLINFITOPTIONS. It has the following fields:
 

int FigureOfMerit;	0: least squares fitting 1: robust fitting, optimizing for minimum absolute deviation (default: 0)
float AbsTolChi;	absolute change of c2 (default: EPSILON)
float FracTolChi;	fractional change of c2 (default: SQRT_EPSILON)
float AbsTolPar;	absolute change of all parameters (default: SQRT_MIN)
float FracTolPar;	fractional change of all parameters (default: SQRT_EPSILON) The four xxxTolChi and xxxTolPar parameters describe the convergence conditions: if the changes achieved in successive iterations are smaller than demanded by these criteria, this signals convergence. Criteria which are not applicable should be set to 0.0.
unsigned HowOftenFulfill;	how often fulfill one of the above conditions before convergence is considered achieved (default: 3)
unsigned LevelOfMethod;	1: Levenberg-Marquardt method, 2: Downhill Simplex (Nelder and Mead) method, 3: both methods alternating; add 4 to this in order to try breaking out of local minima; 0: no fit, calculate only c2 (and Covar) (default: 1)
unsigned LevMarIterations;	max.number of successful iterations of LevMar during one run (default: 100)
unsigned LevMarStarts;	number of LevMar restarts per run (default: 2)
float LambdaStart, LambdaMin, LambdaMax, LambdaDiv, LambdaMul;	treatment of LevMar parameter lambda (don't touch, unless you are an expert!)
unsigned DownhillIterations;	max. number of successful iterations in Downhill Simplex method (default: 200)
float DownhillReflection, DownhillContraction, DownhillExpansion;	treatment of re-shaping of the simplex in Downhill Simplex method (don't touch unless you are an expert!)
unsigned TotalStarts;	max. number of LevMar/Downhill pairs (default: 16)
fVector UpperLimits,         LowerLimits;	impose upper and/or lower limits on parameters. Default for both: NULL or nil
void (*Restrictions)(void);	pointer to a user-defined function, implementing restrictions on the parameters which cannot be formulated as simple upper/lower limits. The function must check the whole parameter vector and edit the parameters as needed. Default: NULL or nil. A typical application of this function would be a model, in which the parameters Ai have to be in a specific order (ascending or descending). So your Restrictions might call VF_sort to enforce that order. Apart from such very special cases, it is generally better not to use Restrictions at all, but let the nonlinfit routine run its course without too much interfering in it.

  A set of options, set by V_setNonlinfitOptions, is valid as long as no other set of options is explicitly passed to VF_nonlinfit. These default options are used by VF_nonlinfit only, if the corresponding parameter FitOptions is set to NULL / nil (or if the simplified syntax is used which does not contain this parameter in the first place).

VF_NONLINFITWORKSPACE

As described above, all non-linear fitting functions need a set of variables for internal use which, for the VF_nonlinfit family of functions, is contained in a struct (Delphi: record) VF_NONLINFITWORKSPACE. It is passed by its address to the respective function. This set of variables need not be initialized. It does not contain user-retrievable information. In the MF_nonlinfit family of functions, a similar set of internal variables is needed as MF_NONLINFITWORKSPACE.
A typical call of VF_nonlinfit would look like this in C/C++:

VF_NONLINFITWORKSPACE ws;
VF_NONLINFITOPTIONS fopt;
V_getNonlinfitOptions( &fopt );
// at this point, modify fopt as desired...
VF_nonlinfit( ParValues, AStatus, nParameters, X, Y, sz, ModelFunc, DerivativeFuncs, &ws, &fopt );

or in Delphi:

ws: VF_NONLINFITWORKSPACE;
fopt: VF_NONLINFITOPTIONS;
V_getNonlinfitOptions( @fopt );
// at this point, modify fopt as desired...
VF_nonlinfit( ParValues, AStatus, nParameters, X, Y, sz, @ModelFunc, @DerivativeFuncs, @ws, @fopt );

If you replace the parameter &ws with NULL / nil, VF_nonlinfit will generate its own workspace.

13.4 Fitting Multiple Data Sets

type VF_EXPERIMENT = record
  X, Y, InvVar: fVector;
  size: UInt;
  WeightOfExperiment: Single;
end;
type PVF_EXPERIMENT = ^VF_EXPERIMENT;
 
type MF_EXPERIMENT = record
  X, Y: fVector;
  MZ, MInvVar: fMatrix;
  htZ, lenZ: UInt;
  WeightOfExperiment: Single;
end;
type PMF_EXPERIMENT = ^MF_EXPERIMENT;

Both in VF_EXPERIMENT and MF_EXPERIMENT, InvVar and WeightOfExperiment are needed only for the weighted variants of the multifit functions.

13.5 Helper Functions for Nonlinear Fits

There are two possibilities how these functions can be called. The first way is to call them from within your model function. For example, you might install a counter into your model function and, upon reaching a certain number of calls, retrieve the current chi2, parameter set, etc. The second way is open only for multithreaded programs: a thread, different from the one performing the fit, may call any of the above functions to control the progress of the fit.
Both of these ways require the VF_NONLINFITWORKSPACE to be globally accessible. When a nonlinfit function is called from multiple threads simultaneously, this poses an obvious problem, especially if you wish to call the helper function from your ModelFunc: How does the model function "know" in real time, from which thread it was called? The solution is: Store a thread index in your parameter vector A. For example, in your main program, you might have something like:

    float A1[11], A2[11], ...;
    int AStatus1[11], AStatus2[11], ...;
    VF_NONLINFITWORKSPACE ws1, ws2, ...;
    ....
    VF_random( A1, 10, 1, -1.0, +1.0 ); // any better guess for starting values would be welcome
    VI_equC( AStatus1, 10, 1 ); // fist 10 parameters are to be fitted
    A1[10] = 1; // thread index = 1
    AStatus1[10] = 0; // of course, this parameter is frozen
    ... throw thread #1 with VF_nonlinfit taking A1, AStatus1, and ws1 as parameters
 
    VF_random( A2, 10, 1, -1.0, +1.0 ); // starting values for second fit
    VI_equC( AStatus2, 10, 1 );
    A2[10] = 2; // thread index = 2
    AStatus2[10] = 0; // again, this parameter is frozen
    ... throw thread #2 with VF_nonlinfit taking A2, AStatus2, and ws2 as parameters

In your model function then, you would have something along the lines of:

void ModelFunc( fVector Y, fVector X, ui size, fVector A )
{
    static unsigned counter1=0, counter2=0, ...;
    float CurrentFOM;
    ....
    switch( (unsigned)(A[10]) )
    {
        case 1: if( (++counter1 % 1024) == 0 )
                {   CurrentFOM = VF_nonlinfit_getFOM( ws1 );
                    printf( "\nCount 1: %u, FOM: %f", counter1, CurrentFOM );
                } break;
        case 1: if( (++counter2 % 1024) == 0 )
                {   CurrentFOM = VF_nonlinfit_getFOM( ws2 );
                    printf( "\nCount 2: %u, FOM: %f", counter2, CurrentFOM );
                } break;
                ....
    }
    ....
}

Back to MatrixLib Table of Contents   OptiVec home

14. Matrix Input and Output

Back to MatrixLib Table of Contents   OptiVec home

15. Graphical Representation of Matrices

Back to MatrixLib Table of Contents   OptiVec home

Copyright © 1996-2022 OptiCode - Dr. Martin Sander Software Development