OptiVec Logo
MatrixLib  



Site Index:

OptiVec home
VectorLib
CMATH
Download
Order
Update
Support

MatrixLib


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
  9. Whole-Matrix Arithmetics: Addition, Multiplication
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.

Back to MatrixLib Table of Contents   OptiVec home

1.1 Matrix Data Types

As VectorLib defines vector data types, here are the matrix data types, defined in MatrixLib:
fMatrixmatrix of floats
dMatrixmatrix of doubles
eMatrixmatrix of extended (long double)
cfMatrixmatrix of fComplex (complex<float>)
cdMatrixmatrix of dComplex (complex<double>)
ceMatrixmatrix of eComplex (complex<extended>)
iMatrixmatrix of int
uMatrixmatrix of unsigned int
and so onfor 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.

Back to MatrixLib Table of Contents   OptiVec home

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
     
    d0u00···  
    l1d1u10··· 
    0l2d2u20···
       ···  
       lN-2dN-2uN-2
       0lN-1dN-1
     
    is compacted into the form
     
    u0 u1 u2 ··· uN-2  *
    d0d1d2···dN-2dN-1 
    * l1l2···lN-2lN-1
     
    The elements l0 and uN-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

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.

Back to MatrixLib Table of Contents   OptiVec home

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.

Back to MatrixLib Table of Contents   OptiVec home


2. Management of Dynamically Generated Matrices

MF_matrixallocate memory for a matrix
MF_matrix0allocate memory and set all elements 0
M_freefree one matrix (data-type independent)
M_nfreefree n matrices (data-type independent; only C/C++)
V_freeAllfree 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_PelementPointer to a specific matrix element
MF_elementvalue 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

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_equ0set all elements to 0
MF_equ1identity matrix: set all diagonal elements to 1.0, all others to 0
MF_equm1negative identity matrix: set all diagonal elements to -1.0, all others to 0
MF_randomfill with random numbers
MF_outerprodmatrix formed by the "outer product" of two vectors
MF_Row_equ0set all elements of one specific row to 0
MF_Col_equ0set all elements of one specific column to 0
MF_Dia_equ0set all diagonal elements to 0
MF_Row_equCset all elements of one specific row to the constant C
MF_Col_equCset all elements of one specific column to the constant C
MF_Dia_equCset all diagonal elements to the constant C
MF_Row_equVcopy a vector into one specific row
MF_Col_equVcopy a vector into one specific column
MF_Dia_equVcopy a vector into the diagonal
MF_Trd_equMcopy a compacted tridiagonal matrix into a general matrix
MF_equMmake one matrix the copy of another
MF_negmake one matrix the negative of another
MCF_conjmake one matrix the complex conjugate of another
MCF_hermconjmake one matrix the Hermitian conjugate of another: MB = MAT*
MF_UequLcopy lower-diagonal elements into upper-diagonal by index-reflection, so as to get a symmetric matrix
MF_LequUcopy upper-diagonal elements into lower-diagonal
 
Two-dimensional windows for spectral analysis are provided by: 
MF_HannHann window
MF_ParzenParzen window
MF_WelchWelch window

Back to MatrixLib Table of Contents   OptiVec home


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_FtoDfMatrix to dMatrix
M_EtoDeMatrix to dMatrix (with overflow protection)
M_CDtoCFcdMatrix to cfMatrix (with overflow protection)
M_DtoEdMatrix to eMatrix

Back to MatrixLib Table of Contents   OptiVec home


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

MF_transposetranspose a matrix: MB = MAT
MCF_hermconjHermitian conjugate: MB = MAT*
MF_rotate90clockwise rotation by 90°
MF_rotate180rotation by 180°
MF_rotate270clockwise rotation by 270° (or counter-clockwise rotation by 90 °)
MF_Rows_revreverse the element order along rows; this corresponds to a reflection of the matrix at the Y axis
MF_Cols_revreverse the element order along columns; this corresponds to a reflection of the matrix at the X axis
MF_Rows_reflectset the upper halves of all rows equal to their reversed lower halves
MF_Cols_reflectset the upper halves of all columns equal to their reversed lower halves
MF_UequLcopy lower-diagonal elements into upper-diagonal by index-reflection, so as to get a symmetric matrix
MF_LequUcopy upper-diagonal elements into lower-diagonal
 
MF_polyinterpolpolynomial interpolation
MF_ratinterpolrational interpolation
MF_natCubSplineInterpolnatural cubic spline interpolation
MF_equMblockextract a block, i.e. a submatrix for which the sampling interval both along rows and along columns is 1
MF_equMblockTextract the transpose of a block
MF_submatrixextract a submatrix with sampling intervals along rows and/or columns possibly different from 1
MF_block_equMcopy a matrix back into a block of another (normally larger) matrix
MF_block_equMTcopy the transpose of a matrix into a block of another (normally larger) matrix
MF_submatrix_equMcopy a submatrix back into another (normally larger) matrix
MF_Row_extractextract a single row and copy it into a vector
MF_Col_extractextract a single column and copy it into a vector
MF_Dia_extractextract the diagonal and copy it into a vector
MF_Trd_extractextract a tridiagonal matrix from a general matrix
MF_Row_insertaugment a matrix by insertion of one row
MF_Col_insertaugment a matrix by insertion of one column
MF_Row_deletedelete one row of a matrix
MF_Col_deletedelete one column of a matrix
 

Back to MatrixLib Table of Contents   OptiVec home


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

MF_Row_negmultiply all elements of a specific row by -1
MF_Col_negmultiply all elements of a specific column by -1
MF_Row_addCadd a constant to all elements of a specific row
MF_Col_addCadd a constant to all elements of a specific column
MF_Dia_addCadd a constant to all diagonal elements
MF_Row_addVadd corresponding vector elements to all elements of a specific row
MF_Col_addVadd corresponding vector elements to all elements of a specific column
MF_Dia_addVadd corresponding vector elements to the diagonal elements
 
A few examples should suffice for the other functions of this family: 
MF_Row_subCsubtract a constant from all elements of a specific row
MF_Col_subrCreverse substraction: difference between column elements and a constant
MF_Dia_mulVmultiply the diagonal elements with corresponding vector elements
MF_Row_divVdivide all elements of a specific row by corresponding vector elements
MF_Col_divrCreverse division: division of a constant by the individual column elements

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

MF_Rows_maxstore the maxima of all rows in a column vector
MF_Cols_maxstore the maxima of all colums in a row vector
MF_Dia_maxreturn the maximum of the diagonal as a scalar
MF_Rows_absmaxstore the absolute maxima of all rows in a column vector
MF_Cols_absmaxstore the absolute maxima of all colums in a row vector
MF_Dia_absmaxreturn the absolute maximum of the diagonal as a scalar
MF_Rows_absminstore the absolute minima of all rows in a column vector
MF_Cols_absminstore the absolute minima of all colums in a row vector
MF_Dia_absminreturn the absolute minimum of the diagonal as a scalar
MF_Rows_sumsum, taken along rows and stored in a column vector
MF_Cols_sumsum, taken along colums and stored in a row vector
MF_Dia_sumsum of the diagonal elements
MF_Rows_prodproduct, taken along rows and stored in a column vector
MF_Cols_prodproduct, taken along colums and stored in a row vector
MF_Dia_prodproduct of the diagonal elements
MF_Rows_runsumrunning sum along rows
MF_Cols_runsumrunning sum along columns
MF_Rows_runprodrunning product along rows
MF_Cols_runprodrunning product along columns
MF_Rows_rotaterotate all rows by a specified number of positions
MF_Cols_rotaterotate all columns by a specified number of positions
MF_Rows_reflectset the upper halves of all rows equal to their reversed lower halves
MF_Cols_reflectset 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_absmaxstore the maxima of the absolute values of all rows in a (real-valued) column vector
MCF_Cols_absmaxstore the maxima of the absolute values of all colums in a (real-valued) row vector
MCF_Dia_absmaxreturn the maximum of the absolute values of the diagonal elements as a (real) scalar
MCF_Rows_absmaxReImfind 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_absmaxReImfind 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_absmaxReImfind 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_cabsmaxfind the complex numbers of largest magnitude along rows and store these row maxima in a (complex) column vector
MCF_Cols_cabsmaxfind the complex numbers of largest magnitude along columns and store these column maxima in a (complex) row vector
MCF_Dia_cabsmaxfind the complex number of largest magnitude along the diagonal of a square matrix and return it as a (complex) scalar
MCF_Rows_maxReImfind 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_maxReImfind 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_maxReImfind 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_sabsmaxfind the complex numbers of largest sum |Re|+|Im| along rows and store these row maxima in a (complex) column vector
MCF_Cols_sabsmaxfind the complex numbers of largest sum |Re|+|Im| along columns and store these column maxima in a (complex) row vector
MCF_Dia_sabsmaxfind 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_centerOfGravityIndcenter of gravity, returned as interpolated element indices
MF_centerOfGravityVcenter 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

MF_Rows_exchangeexchange two rows
MF_Cols_exchangeexchange two columns
MF_Rows_addadd one row to another (destination += source)
MF_Cols_addadd one column to another
MF_Rows_subsubtract one row from another (destination -= source)
MF_Cols_subsubtract one column from another
MF_Rows_Caddadd scaled row to another (destination += C * source)
MF_Cols_Caddadd scaled column to another
MF_Rows_lincomblinear combination of two rows
MF_Cols_lincomblinear combination of two columns

Back to MatrixLib Table of Contents   OptiVec home


9. Whole-Matrix Arithmetics: Addition, Multiplication

a) Element-wise operations
MF_addMadd two matrices
MF_addMTadd one matrix and the transpose of another matrix
MC = MA + MBT
MF_subMsubtract one matrix from another
MF_subMTsubtract a transposed matrix
MC = MA - MBT
MF_subrMTsubtract a matrix from another, transposed, matrix
MC = MBT - MA
MF_mulCmultiply all matrix elements by a constant
MCF_mulReCmultiply all elements of a complex matrix by a real number
MF_divCdivide all matrix elements by a constant
MCF_divReCdivide all elements of a complex matrix by a real number
MFs_addMscaled addition of two matrices:
MC = c * (MA + MB)
MFs_addMTscaled addition of one matrix and the transpose of another:
MC = c * (MA + MBT)
MFs_subMscaled subtraction of two matrices:
MC = c * (MA - MB)
MFs_subMTscaled subtraction of one matrix and the transpose of another:
MC = c * (MA - MBT)
MFs_subrMTscaled reverse subtraction of one matrix and the transpose of another:
MC = c * (MBT - MA)
MF_lincomblinear combination:
MC = ca * MA + cb * MB

b)Matrix multiplication:
MF_mulVmultiply a matrix by a column vector:
Y = MA * X
MF_TmulVmultiply the transpose of a matrix by a column vector:
Y = MAT * X
VF_mulMmultiply a row vector by a matrix:
Y = X * MA
VF_mulMTmultiply a row vector by the transpose of a matrix:
Y = X * MAT
MF_mulMmultiply two matrices:
MC = MA * MB
MF_mulMTmultiply one matrix by the transpose of another matrix:
MC = MA * MBT
MF_TmulMmultiply the transpose of a matrix by another matrix:
MC = MAT * MB
MF_TmulMTmultiply the transpose of one matrix by the transpose of another matrix:
MC = MAT * MBT
MCF_mulMHmultiply one matrix by the hermitian conjugate of another matrix:
MC = MA * MBT *
MCF_HmulMmultiply 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_mulMdiamultiply a general matrix by a diagonal matrix:
MC = MA * MBDia
MF_TmulMdiamultiply the transpose of a matrix by a diagonal matrix:
MC = MAT * MBDia
MFdia_mulMmultiply a diagonal matrix by a general matrix:
MC = MADia * MB
MFdia_mulMTmultiply a diagonal matrix by the transpose of a general matrix:
MC = MADia * MBT

Back to MatrixLib Table of Contents   OptiVec home


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

Back to MatrixLib Table of Contents   OptiVec home


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

Back to MatrixLib Table of Contents   OptiVec home


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_FFTtoCForward Fast Fourier Transform (FFT) of a real matrix; the result is a cartesian complex matrix
MF_FFTForward 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_FFTForward and backward FFT of complex matrices
MF_convolveConvolution with a spatial response function
MF_deconvolveDeconvolution
MF_filterSpatial filtering
MF_autocorrSpatial autocorrelation
MF_xcorrSpatial cross-correlation
MF_spectrumSpatial frequency spectrum

The one-dimensional Fourier Transform along rows or along columns is available as:
MF_Rows_FFTFourier Transform along rows
MF_Cols_FFTFourier Transform along columns

Back to MatrixLib Table of Contents   OptiVec home


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_linfitfit X-Y-Z data to a function linear in its parameters
MF_linfitwWthe same with non-equal weighting of individual data points
MF_nonlinfitfit X-Y-Z data to an arbitrary, possibly non-linear function
MF_nonlinfitwWthe same for non-equal data-point weighting
MF_multiLinfitfit multiple X-Y-Z data sets to one common linear function
MF_multiLinfitwWthe same for non-equal data-point weighting
MF_multiNonlinfitfit 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_getBestValuesbest set of parameters found so far
VF_nonlinfit_getFOMthe best figure-of-merit (c2, chi-square, or |c|, chiabs, depending on the VF_NONLINFITOPTIONS member "FigureOfMerit") obtained so far
VF_nonlinfit_getTestDirreturns 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_getTestParindex of the parameter currently under "breakout" investigation
VF_nonlinfit_getTestRunindex of the current "breakout" test run (for each fitted parameter, one test run is performed)
VF_nonlinfit_stopstops 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;
                ....
    }
    ....
}

Back to MatrixLib Table of Contents   OptiVec home

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_printprint a matrix to the screen (without paging or row cut-off); (DOS, EasyWin and Delphi only)
MF_fprintprint a matrix in ASCII format to a stream
MF_storestore in binary format
MF_recallretrieve in binary format
MF_writewrite in ASCII format in a stream
MF_readread from an ASCII file

Back to MatrixLib Table of Contents   OptiVec home


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_zAutoDensityMapColor 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_zDataDensityMapColor 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_setDensityBoundsSet a color scale for matrix color-density plots.
M_setDensityMapBoundsSet a color scale and draw an X-Y coordinate system for matrix color-density plots.
M_findDensityMapBoundsCalculate 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.

Back to MatrixLib Table of Contents   OptiVec home

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

Last modified: 21 June 2022