OptiVec home
VectorLib
CMATH
Order
Update
Support

# Contents

1. Introduction
1.1 Matrix Data Types
1.2 Prefixes of MatrixLib Functions
1.3 C/C++ Specifics
1.3.1 MatObj, the object-oriented interface for matrix functions
1.3.2 Direct-pointer form of matrix function calls
1.4 Pascal/Delphi Specifics
2. Management of Dynamically Generated Matrices
3. Initialization of Matrices
4. Data-Type Conversions
5. Symmetry Operations (Transposing, Rotating, Reflecting),  Interpolation, Augmenting, Deleting, Extracting, and Filling Parts of a Matrix
6. Arithmetic Operations Performed on a Single Row, Column, or the Diagonal
7. Operations Performed Along All Rows or All Columns Simultaneously, or the Diagonal; Center of Gravity of a Matrix
8. Operations Involving Two Rows or Two Colums
10. Linear Algebra
11. Eigenvalues and Eigenvectors, Matrix Square-Root
12. Two-Dimensional Fourier-Transform Methods
13. Data Fitting
13.1 Polynomials
13.2 General Linear Model Functions
13.3 Non-Linear Models
13.4 Fitting Multiple Data Sets
13.5 Helper Functions for Nonlinear Fits
13.6 Summary of Fitting Functions
14. Matrix Input and Output
15. Graphical Representation of Matrices

## 1. Introduction

The MatrixLib part of OptiVec builds upon the principles established in the VectorLib part. You may wish to refer to the introductory chapters of the VectorLib documentation, before reading on for MatrixLib.

### 1.1 Matrix Data Types

As VectorLib defines vector data types, here are the matrix data types, defined in MatrixLib:
 fMatrix matrix of floats dMatrix matrix of doubles eMatrix matrix of extended (long double) cfMatrix matrix of fComplex (complex) cdMatrix matrix of dComplex (complex) ceMatrix matrix of eComplex (complex) iMatrix matrix of int uMatrix matrix of unsigned int and so on for all other integer 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.

### 1.2 Prefixes of MatrixLib Functions

Each MatrixLib function has a prefix defining the data type on which it acts:
 Prefix: Arguments: MF_ fMatrix, float and fVector MD_ dMatrix, double and dVector ME_ eMatrix, extended (long double) and eVector MCF_ cfMatrix, fComplex and cfVector MCD_ cdMatrix, dComplex and cdVector MCE_ ceMatrix, eComplex and ceVector MI_ iMatrix, int and iVector MU_ uMatrix, unsigned int and uVector similarly for all other integer types

In a few cases, the prefix is augmented by a three-letter code denoting special matrix properties:
• 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 l0 and uN-1, marked with an asterisk, are undefined and never used.

### 1.3.1 MatObj, the object-oriented interface for matrix functions

Similarly to the vector objects of VecObj, the object-oriented interface for the matrix functions is contained in the include-files <fMatObj.h>, <dMatObj.h> etc., with one include-file for each of the data-types supported in OptiVec. Remember that, in order to get the whole interface (for all data types at once),
#include <OptiVec.h>.
For access to any of the vector graphics functions, always include <OptiVec.h>.

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

In the C/C++ version of MatrixLib, all functions taking matrix arguments exist in a second form. In this form, all matrix arguments are replaced by pointers to the first elements of the respective matrices. You can explicitly use the second form by leaving away the underbar of the function-name prefix. Calling
MF_mulM( MC, MA, MB, htA, lenA, lenB );
is equivalent to calling
MFmulM(&(MC[0][0]),&(MA[0][0]),&(MB[0][0]),htA,lenA,lenB); or
MFmulM( MC[0], MA[0], MB[0], htA, lenA, lenB );
Actually, the run-time routines are in the second form, and the macros defined in <MFstd.h> etc. convert calls to the first form into calls to the second form.

### 1.4 Pascal/Delphi Version Specifics:

In the Pascal/Delphi version of MatrixLib, matrices of all data types are actually defined as pointers to the element M[0][0]. This means you can pass static matrices (like MX: array[0..5][0..5] of Single; ) to MatrixLib functions with the address operator:
MF_equ1( @(MX[0][0]), 6 );

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.

## 2. Management of Dynamically Generated Matrices

 MF_matrix allocate memory for a matrix MF_matrix0 allocate memory and set all elements 0 M_free free one matrix (data-type independent) M_nfree free n matrices (data-type independent; only C/C++) V_freeAll free all existing vectors and matrices

OptiVec's dynamically allocated matrices are aligned on 32-byte boundaries, which allows for optimum cache-line matching for the Pentium processor and its currently available successors and competitors.

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 *)

## 3. Initialization of Matrices

In order to initialize all matrix elements with the same value, or to perform the same arithmetic operation on all matrix elements simultaneously, please use the respective VectorLib function. You can do so, because all matrices are in fact stored as vectors, occupying a contiguous space in memory. All you have to do is to pass the first row of the matrix (rather than the matrix itself) to the vector function. Fore example, you can initialize all matrix elements with a constant C by calling
VF_equC( MA[0], ht*len, C ); /* C/C++ */
VF_equC( MA, ht*len, C ); (* Pascal/Delphi *)

As you shall see, some of the most common operations of this kind are also explicitly defined in MatrixLib, like initialization with zero, MF_equ0, or multiplication by a constant, available as MF_mulC (see chapter 9). Here are the typical matrix initialization functions:

 MF_equ0 set all elements to 0 MF_equ1 identity matrix: set all diagonal elements to 1.0, all others to 0 MF_equm1 negative identity matrix: set all diagonal elements to -1.0, all others to 0 MF_random fill with random numbers MF_outerprod matrix formed by the "outer product" of two vectors MF_Row_equ0 set all elements of one specific row to 0 MF_Col_equ0 set all elements of one specific column to 0 MF_Dia_equ0 set all diagonal elements to 0 MF_Row_equC set all elements of one specific row to the constant C MF_Col_equC set all elements of one specific column to the constant C MF_Dia_equC set all diagonal elements to the constant C MF_Row_equV copy a vector into one specific row MF_Col_equV copy a vector into one specific column MF_Dia_equV copy a vector into the diagonal MF_Trd_equM copy a compacted tridiagonal matrix into a general matrix MF_equM make one matrix the copy of another MF_neg make one matrix the negative of another MCF_conj make one matrix the complex conjugate of another MCF_hermconj make one matrix the Hermitian conjugate of another: MB = MAT* MF_UequL copy lower-diagonal elements into upper-diagonal by index-reflection, so as to get a symmetric matrix MF_LequU copy upper-diagonal elements into lower-diagonal

Two-dimensional windows for spectral analysis are provided by:
 MF_Hann Hann window MF_Parzen Parzen window MF_Welch Welch window

## 4. Data-Type Conversions

Matrices of every data type can be converted into every other. Only a few examples are given; the rest should be obvious.
 M_FtoD fMatrix to dMatrix M_EtoD eMatrix to dMatrix (with overflow protection) M_CDtoCF cdMatrix to cfMatrix (with overflow protection) M_DtoE dMatrix to eMatrix

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

 MF_transpose transpose a matrix: MB = MAT MCF_hermconj Hermitian conjugate: MB = MAT* MF_rotate90 clockwise rotation by 90° MF_rotate180 rotation by 180° MF_rotate270 clockwise rotation by 270° (or counter-clockwise rotation by 90 °) MF_Rows_rev reverse the element order along rows; this corresponds to a reflection of the matrix at the Y axis MF_Cols_rev reverse the element order along columns; this corresponds to a reflection of the matrix at the X axis MF_Rows_reflect set the upper halves of all rows equal to their reversed lower halves MF_Cols_reflect set the upper halves of all columns equal to their reversed lower halves MF_UequL copy lower-diagonal elements into upper-diagonal by index-reflection, so as to get a symmetric matrix MF_LequU copy upper-diagonal elements into lower-diagonal

 MF_polyinterpol polynomial interpolation MF_ratinterpol rational interpolation MF_natCubSplineInterpol natural cubic spline interpolation MF_equMblock extract a block, i.e. a submatrix for which the sampling interval both along rows and along columns is 1 MF_equMblockT extract the transpose of a block MF_submatrix extract a submatrix with sampling intervals along rows and/or columns possibly different from 1 MF_block_equM copy a matrix back into a block of another (normally larger) matrix MF_block_equMT copy the transpose of a matrix into a block of another (normally larger) matrix MF_submatrix_equM copy a submatrix back into another (normally larger) matrix MF_Row_extract extract a single row and copy it into a vector MF_Col_extract extract a single column and copy it into a vector MF_Dia_extract extract the diagonal and copy it into a vector MF_Trd_extract extract a tridiagonal matrix from a general matrix MF_Row_insert augment a matrix by insertion of one row MF_Col_insert augment a matrix by insertion of one column MF_Row_delete delete one row of a matrix MF_Col_delete delete one column of a matrix

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

A few examples should suffice for the other functions of this family:
 MF_Row_subC subtract a constant from all elements of a specific row MF_Col_subrC reverse substraction: difference between column elements and a constant MF_Dia_mulV multiply the diagonal elements with corresponding vector elements MF_Row_divV divide all elements of a specific row by corresponding vector elements MF_Col_divrC reverse division: division of a constant by the individual column elements

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

 MF_Rows_max store the maxima of all rows in a column vector MF_Cols_max store the maxima of all colums in a row vector MF_Dia_max return the maximum of the diagonal as a scalar MF_Rows_absmax store the absolute maxima of all rows in a column vector MF_Cols_absmax store the absolute maxima of all colums in a row vector MF_Dia_absmax return the absolute maximum of the diagonal as a scalar MF_Rows_absmin store the absolute minima of all rows in a column vector MF_Cols_absmin store the absolute minima of all colums in a row vector MF_Dia_absmin return the absolute minimum of the diagonal as a scalar MF_Rows_sum sum, taken along rows and stored in a column vector MF_Cols_sum sum, taken along colums and stored in a row vector MF_Dia_sum sum of the diagonal elements MF_Rows_prod product, taken along rows and stored in a column vector MF_Cols_prod product, taken along colums and stored in a row vector MF_Dia_prod product of the diagonal elements MF_Rows_runsum running sum along rows MF_Cols_runsum running sum along columns MF_Rows_runprod running product along rows MF_Cols_runprod running product along columns MF_Rows_rotate rotate all rows by a specified number of positions MF_Cols_rotate rotate all columns by a specified number of positions MF_Rows_reflect set the upper halves of all rows equal to their reversed lower halves MF_Cols_reflect set the upper halves of all columns equal to their reversed lower halves

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

## 8. Operations Involving Two Rows or Two Colums

 MF_Rows_exchange exchange two rows MF_Cols_exchange exchange two columns MF_Rows_add add one row to another (destination += source) MF_Cols_add add one column to another MF_Rows_sub subtract one row from another (destination -= source) MF_Cols_sub subtract one column from another MF_Rows_Cadd add scaled row to another (destination += C * source) MF_Cols_Cadd add scaled column to another MF_Rows_lincomb linear combination of two rows MF_Cols_lincomb linear combination of two columns

## 9. Whole-Matrix Arithmetics: Addition, Multiplication

a) Element-wise operations
 MF_addM add two matrices MF_addMT add one matrix and the transpose of another matrix MC = MA + MBT MF_subM subtract one matrix from another MF_subMT subtract a transposed matrix MC = MA - MBT MF_subrMT subtract a matrix from another, transposed, matrix MC = MBT - MA MF_mulC multiply all matrix elements by a constant MCF_mulReC multiply all elements of a complex matrix by a real number MF_divC divide all matrix elements by a constant MCF_divReC divide all elements of a complex matrix by a real number MFs_addM scaled addition of two matrices: MC = c * (MA + MB) MFs_addMT scaled addition of one matrix and the transpose of another:MC = c * (MA + MBT) MFs_subM scaled subtraction of two matrices: MC = c * (MA - MB) MFs_subMT scaled subtraction of one matrix and the transpose of another:MC = c * (MA - MBT) MFs_subrMT scaled reverse subtraction of one matrix and the transpose of another:MC = c * (MBT - MA) MF_lincomb linear combination: MC = ca * MA + cb * MB

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

## 10. Linear Algebra

There are four groups of linear algebra functions. The first group consists of "easy-to-use" versions which can be used as black-box functions without caring about their internal working. The second group consists of functions for LU decomposition and its applications. The third group consists of Cholesky factorization and its applications on symmetric, positive-definite matrices. The fourth group, finally, is devoted to Singular Value Decomposition (SVD). The functions based on LUD are available both for real and for complex matrices, those based on Cholesky and SVD only for real matrices.
Here are the "easy-to-use" versions of linear algebra functions:
 solve simultaneous linear equations (using LU decomposition) invert a matrix determinant of a matrix solve a symetric system of linear equation, assuming the system is positive-definite (using Cholesky decomposition; in case the system turns out not to be positive-definite, LU decomposition is used) invert a symmetric matrix (trying first Cholesky decomposition, but switching to LUD in case the system turns out not to be positive-definite) determinant of a symmetric matrix solve simultaneous linear equations, using Singular Value Decomposition tries first solution by LUD; if that fails, SVD is done

Now some functions for explicit LU decomposition and for treatment of LU decomposed matrices:
 decompose into LU form check if was successful set editing threshold for ; may be used to work around singularities retrieve currently set threshold solve simultaneous linear equations, given the matrix in LU form improve the accuracy of the solution of an LU-decomposed linear system by iteration invert matrix already composed into LU form determinant of matrix already composed into LU form

The following functions are available for Cholesky decomposition and its applications on positive-definite matrices:
 Cholesky decomposition into left-triangular form Cholesky decomposition into right-triangular form solution of a positive-definite system of linear equations after Cholesky decomposition into left-triangular form solution of a positive-definite system of linear equations after Cholesky decomposition into right-triangular form iterative improvement of the accuracy of a solution obtained via Cholesky decomposition into L form iterative improvement of the accuracy of a solution obtained via Cholesky decomposition into R form invert a matrix which has already been reduced into L form by Cholesky decomposition invert a matrix which has already been reduced into L form by Cholesky decomposition determinant of a matrix which has already been reduced into L or R form by Cholesky decomposition

Singular Value Decomposition and related functions are offered as:
 Singular Value Decomposition Sorting of Singular Values in descending order with comcomitant re-ordering of left and right singular vectors solve SV-decomposed set of linear equations set threshold for Singular Value editing retrieve current threshold for Singular Value editing

## 11. Eigenvalues and Eigenvectors, Matrix Square-Root

At present, only the special, but most frequent case of a symmetric real matrix is covered:
 MFsym_eigenvalues eigenvalues with or without eigenvectors of a symmetric real matrix MFsym_sqrt square-root of a symmetric, positive definite matrix

## 12. Two-Dimensional Fourier-Transform Methods

By analogy with the corresponding one-dimensional Fourier Transform methods described in chapter 4.8 of the VectorLib documentation, the following functions for the two-dimensional case are available in MatrixLib:
 MF_FFTtoC Forward Fast Fourier Transform (FFT) of a real matrix; the result is a cartesian complex matrix MF_FFT Forward and backward FFT of real matrices; using symmetry relations, the complex result is packed into a real matrix of the same size as the input matrix MCF_FFT Forward and backward FFT of complex matrices MF_convolve Convolution with a spatial response function MF_deconvolve Deconvolution MF_filter Spatial filtering MF_autocorr Spatial autocorrelation MF_xcorr Spatial cross-correlation MF_spectrum Spatial frequency spectrum

The one-dimensional Fourier Transform along rows or along columns is available as:
 MF_Rows_FFT Fourier Transform along rows MF_Cols_FFT Fourier Transform along columns

## 13. Data Fitting

MatrixLib provides a broad range of data fitting functions for various classes of model functions. Let us first show you an overview over the available functions, before we describe the various classes of model functions and their treatment in OptiVec fitting functions in detail.

 VF_linregress equally-weighted linear regression on X-Y data VF_linregresswW the same with non-equal weighting VF_polyfit fitting of one X-Y data set to a polynomial VF_polyfitwW the same for non-equal data-point weighting VF_linfit fitting of one X-Y data set to an arbitrary function linear in its parameters VF_linfitwW the same for non-equal data-point weighting VF_nonlinfit fitting of one X-Y data set to an arbitrary, possibly non-linear function VF_nonlinfitwW the same for non-equal data-point weighting VF_multiLinfit fitting of multiple X-Y data sets to one common linear function VF_multiLinfitwW the same for non-equal data-point weighting VF_multiNonlinfit fitting of multiple X-Y data sets to one common nonlinear function VF_multiNonlinfitwW the same for non-equal data-point weighting MF_linfit fit X-Y-Z data to a function linear in its parameters MF_linfitwW the same with non-equal weighting of individual data points MF_nonlinfit fit X-Y-Z data to an arbitrary, possibly non-linear function MF_nonlinfitwW the same for non-equal data-point weighting MF_multiLinfit fit multiple X-Y-Z data sets to one common linear function MF_multiLinfitwW the same for non-equal data-point weighting MF_multiNonlinfit fit multiple X-Y-Z data sets to one common nonlinear function MF_multiNonlinfitwW the same for non-equal data-point weighting

Now let us have a look at the various classes of model functions. The most simple one is a straight line,
Yi = a * Xi + b;
Fitting Y = f(X) data to a straight line is called "linear regression" and accomplished by VF_linregress.
The next level of sophistication are polynomials, described in the following paragraph:

### 13.1 Polynomials

Polynomials are functions of the form
Yi = a0 + a1Xi + a2Xi2 ... anXin
Just as linear regression, polynomial fitting is available for Y-X data (VF_polyfit), but not for MZ=f(X,Y) data. In polynomial fitting, you are interested in the coefficients ai up to a certain limit, the degree n of the polynomial. In the simplest form, all coefficients are treated as free parameters, without the possibility of "fixing" those whose values you happen to know beforehand. Polynomial fitting is most useful for degrees of 2 to 4. With degrees higher than 5, you will be able to fit almost any data - but the resulting coefficients will be at best inaccurate and at worst useless, as already very little experimental noise will strongly influence the balance of the higher terms. If you do need polynomials of higher degrees, you should carefully examine your problem to see which terms you can eliminate. This leads us from polynomials to the next class of fitting functions, namely:

### 13.2 General Linear Model Functions

In the general linear case, you have to write the model function yourself and pass it as an argument to the fitting routine. The word "linear" means that the model function consists of a linear combination of basis functions each of which depends linearly on the fitting coefficients:
y = a0f0(x) + a1f1(x) + a2f2(x)...
The individual functions fi(x) may be non-linear in x, but not in the coefficients ai. For example,
y = a0sin(x)+ a1cos(x)
is a valid linear fitting function, but
y = sin(a0x)+ cos(a1x)
is not.
For general linear fits, the model function has to be provided in a form which calculates the individual basis functions for each x-value present in the data set. The above example with sine and cosine functions would be coded in C/C++ as

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

Note that the coefficients ai 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 statusi is set to 0, the corresponding parameter, ai, 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 ai 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(a0x)+ cos(a1x)
require their own treatment, we arrive at the last and most general class of model functions:

### 13.3 Non-Linear Models

As described above, fits to linear model functions (which include, of course, simple polynomials) are performed internally by solving a set of coupled linear equations - which is basically a simple matrix inversion process. The situation for non-linear model functions is completely different: no closed-form solution exists, and the fitting problem can be solved only iteratively. Therefore, fitting to non-linear models is much more time-consuming (by 3-5 orders of magnitude!) than fitting to linear models. Two non-linear fitting algorithms are offered:
• the Levenberg-Marquardt method, and the
• Downhill-Simplex method by Nelder and Mead.
As a starting-point, it is normally a good idea to choose a combination of both (choose FitOptions.LevelOfMethod = 3, see below).
The fundamental difference between linear and non-linear data fitting is reflected also in the required formulation of the model functions. In contrast to the linear case, here it is not the individual basis functions which are needed. Rather, the model function will be called by VF_nonlinfit with the whole X-vector as input argument (as well as its size and the array A of parameters) and has to return the whole Y-vector. For the above non-linear sine and cosine example, this could be coded in C/C++ as

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 ai, you could also calculate dY / dai 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

In addition to the fitting functions for single data sets, there are routines for fitting multiple data sets simultaneously. Say, you are a physicist or a chemist and want to measure a rate constant k. You have performed the same kinetic measurement ten times under slightly different conditions. All these measurements have, of course, the same k, but other parameters, like amplitude and time-zero, will differ from experiment to experiment. Now, the usual way would be to perform one fit for each of these ten data sets and average over the resulting k's.
There is a problem with this approach: you don't really know which weight you have to assign to each of the measurements, and how to treat the overlapping error margins of the individual fits. You might want to perform a fit on all your data sets simultaneously, taking the same k for all data sets, and assigning individual amplitudes etc. to each set. This is precisely what the multiLinfit and multiNonlinfitfunctions of MatrixLib are designed for. You can fit multiple data-sets both to linear and nonlinear models, for Y = f(X) data as well as for MZ = f(X,Y) data.
The multiLinfit and multiNonlinfit functions need the input data to be passed in structures named VF_EXPERIMENT for X-Y data and MF_EXPERIMENT for X-Y-Z data. In C/C++, VF_EXPERIMENT has the following fields:

 fVector X, Y; X and Y vectors (independent variables) fMatrix MZ; measured data z=f(x,y) fMatrix MInvVar; inverse variances of the individual matrix elements (needed only for the "with weights" functions) unsigned htZ, lenZ; matrix dimensions float WeightOfExperiment; individual weight to be assigned to the whole experiment (again needed only for the weighted variants)

In C/C++, MF_EXPERIMENT has the following fields:

 fVector X, Y; X and Y vectors fVector InvVar; inverse variances of the individual vector elements (needed only for the "with weights" functions) ui size; the number of vector elements float WeightOfExperiment; individual weight to be assigned to the whole experiment (again needed only for the weighted variants)

In Pascal/Delphi, VF_EXPERIMENT and MF_EXPERIMENT are defined as follows:

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

As nonlinear fitting tasks - especially with multiple data sets - may become quite tedious and take a very long time (sometimes many hours!) to converge, there is a number of functions which allow to follow the course of nonlinear fits. The names of these functions are always derived from the fitting function they are used with. For example, the helper functions for VF_nonlinfit are:

 VF_nonlinfit_getBestValues best set of parameters found so far VF_nonlinfit_getFOM the best figure-of-merit (c2, chi-square, or |c|, chiabs, depending on the VF_NONLINFITOPTIONS member "FigureOfMerit") obtained so far VF_nonlinfit_getTestDir returns the test direction (+1 for upwards, -1 for downwards) during "breakout" attempts (level-of-method greater than 3, see the description of VF_NONLINFITOPTIONS below) VF_nonlinfit_getTestPar index of the parameter currently under "breakout" investigation VF_nonlinfit_getTestRun index of the current "breakout" test run (for each fitted parameter, one test run is performed) VF_nonlinfit_stop stops fitting after completion of the current Levenberg-Marquardt or Downhill cycle

For the other nonlinear fitting functions, the corresponding helper-function names are obtained by replacing the "VF_nonlinfit_" prefix with the respective fitting function name, as in VF_nonlinfitwW_getBestValues, MF_multiNonlinfit_stop, and so on.

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;
....
}
....
}

## 14. Matrix Input and Output

The matrix input/output functions are all analogous to the corresponding vector functions
 MF_cprint print a matrix to the screen. If necessary, rows are cut off at the screen boundaries. If there are more rows than screen lines, proper paging is applied. (DOS and Delphi only) MF_print print a matrix to the screen (without paging or row cut-off); (DOS, EasyWin and Delphi only) MF_fprint print a matrix in ASCII format to a stream MF_store store in binary format MF_recall retrieve in binary format MF_write write in ASCII format in a stream MF_read read from an ASCII file

## 15. Graphical Representation of Matrices

True 3D-plotting functions will be included in future versions. For now, only color-density plots are available. In these plots, each data value is translated into a color value by linear interpolation between two colors, specified as mincolor and maxcolor.

 MF_xyzAutoDensityMap Color density map for z=f(x,y) with automatic scaling of the x and y axes and of the color density scale between mincolor and maxcolor MF_xyzDataDensityMap z=f(x,y) color density map, plotted into an existing axis frame, and using the color density scale set by the last call to an AutoDensityMap function MF_zAutoDensityMap Color density map for z=f(i,j) with automatic scaling of the x and y axes and of the color density scale between mincolor and maxcolor. i and j are the indices in x and y direction, respectively MF_zDataDensityMap Color density map for z=f(i,j), plotted into an existing axis frame, and using the color density scale set by the last call to an AutoDensityMap function M_setDensityBounds Set a color scale for matrix color-density plots. M_setDensityMapBounds Set a color scale and draw an X-Y coordinate system for matrix color-density plots. M_findDensityMapBounds Calculate a color scale and draw an X-Y coordinate system for matrix color-density plots, ensuring that the grid lines of the coordinate system correspond to exact (rather than only rounded) values.