
MatrixLibContents1. Introduction1.1 Matrix Data Types 1.2 Prefixes of MatrixLib Functions 1.3 C/C++ Specifics 1.3.1 MatObj, the objectoriented interface for matrix functions 1.3.2 Directpointer form of matrix function calls 1.4 Pascal/Delphi Specifics 2. Management of Dynamically Generated Matrices 3. Initialization of Matrices 4. DataType 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. WholeMatrix Arithmetics: Addition, Multiplication 10. Linear Algebra 11. Eigenvalues and Eigenvectors, Matrix SquareRoot 12. TwoDimensional FourierTransform Methods 13. Data Fitting 13.1 Polynomials 13.2 General Linear Model Functions 13.3 NonLinear 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
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.

fMatrix  matrix of floats 
dMatrix  matrix of doubles 
eMatrix  matrix of extended (long double) 
cfMatrix  matrix of fComplex (complex<float>) 
cdMatrix  matrix of dComplex (complex<double>) 
ceMatrix  matrix of eComplex (complex<extended>) 
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 twodimensional arrays provided by the respective target compilers. This means that the matrices are stored rowwise in MatrixLib versions for C and C++ compilers, but columnwise 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
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 
d_{0}  u_{0}  0  ···  
l_{1}  d_{1}  u_{1}  0  ···  
0  l_{2}  d_{2}  u_{2}  0  ··· 
···  
l_{N2}  d_{N2}  u_{N2}  
0  l_{N1}  d_{N1} 
u_{0}  u_{1}  u_{2}  ···  u_{N2}  * 
d_{0}  d_{1}  d_{2}  ···  d_{N2}  d_{N1} 
*  l_{1}  l_{2}  ···  l_{N2}  l_{N1} 
Back to MatrixLib Table of Contents OptiVec home
By analogy with the alias names of the vector objects, the matrix objects get the alias names fMatObj, dMatObj, and so on, with the datatype 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 vectorvector operator * refers to elementwise multiplication, the matrixmatrix, vectormatrix, or matrixvector operator * means true matrixmatrix or matrixvector multiplication.
If you ever need to process a MatObj matrix in a "classic" plainC 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:
Back to MatrixLib Table of Contents OptiVec home
The situation is somewhat different for dynamic arrays in Delphi. While, in the onedimensional case, they they can be used with all VectorLib functions by simply typecasting fArrays into fVectors, twodimensional arrays require a different approach: Dynamic twodimensional arrays can be created in Delphi by declaring them as "array of fArray", "array of dArray", etc. Shortcut definitions for these types are also given in the VecLib unit, as f2DArray, d2DArray, etc. As a consequence of the described way of defining 2DArrays 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
MF_matrix  allocate memory for a matrix 
MF_matrix0  allocate memory and set all elements 0 
M_free  free one matrix (datatype independent) 
M_nfree  free n matrices (datatype independent; only C/C++) 
V_freeAll  free all existing vectors and matrices 
C/C++ version only:
OptiVec's dynamically allocated matrices can be addressed just like twodimensional 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 
Back to MatrixLib Table of Contents OptiVec home
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 = MA^{T*} 
MF_UequL  copy lowerdiagonal elements into upperdiagonal by indexreflection, so as to get a symmetric matrix 
MF_LequU  copy upperdiagonal elements into lowerdiagonal 
MF_Hann  Hann window 
MF_Parzen  Parzen window 
MF_Welch  Welch window 
Back to MatrixLib Table of Contents OptiVec home
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 
Back to MatrixLib Table of Contents OptiVec home
MF_transpose  transpose a matrix: MB = MA^{T} 
MCF_hermconj  Hermitian conjugate: MB = MA^{T*} 
MF_rotate90  clockwise rotation by 90° 
MF_rotate180  rotation by 180° 
MF_rotate270  clockwise rotation by 270° (or counterclockwise 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 lowerdiagonal elements into upperdiagonal by indexreflection, so as to get a symmetric matrix 
MF_LequU  copy upperdiagonal elements into lowerdiagonal 
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 
Back to MatrixLib Table of Contents OptiVec home
MF_Row_neg  multiply all elements of a specific row by 1 
MF_Col_neg  multiply all elements of a specific column by 1 
MF_Row_addC  add a constant to all elements of a specific row 
MF_Col_addC  add a constant to all elements of a specific column 
MF_Dia_addC  add a constant to all diagonal elements 
MF_Row_addV  add corresponding vector elements to all elements of a specific row 
MF_Col_addV  add corresponding vector elements to all elements of a specific column 
MF_Dia_addV  add corresponding vector elements to the diagonal elements 
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 
Back to MatrixLib Table of Contents OptiVec home
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 (realvalued) column vector 
MCF_Cols_absmax  store the maxima of the absolute values of all colums in a (realvalued) row vector 
MCF_Dia_absmax  return the maximum of the absolute values of the diagonal elements as a (real) scalar 
MCF_Rows_absmaxReIm  find the maximum absolute values of all real and of all imaginary parts separately along the rows of a matrix; merge the maxima into complex numbers, and store them in a column vector 
MCF_Cols_absmaxReIm  find the maximum absolute values of all real and of all imaginary parts separately along the columns of a matrix; merge the maxima into complex numbers, and store them in a row vector 
MCF_Dia_absmaxReIm  find the maximum absolute values of all real and of all imaginary parts separately along the diagonal of a square matrix; merge these two maxima into one complex number and return it as a scalar 
MCF_Rows_cabsmax  find the complex numbers of largest magnitude along rows and store these row maxima in a (complex) column vector 
MCF_Cols_cabsmax  find the complex numbers of largest magnitude along columns and store these column maxima in a (complex) row vector 
MCF_Dia_cabsmax  find the complex number of largest magnitude along the diagonal of a square matrix and return it as a (complex) scalar 
MCF_Rows_maxReIm  find the maxima of all real and of all imaginary parts separately along the rows of a matrix; merge the maxima into complex numbers, and store them in a column vector 
MCF_Cols_maxReIm  find the maxima of all real and of all imaginary parts separately along the columns of a matrix; merge the maxima into complex numbers, and store them in a row vector 
MCF_Dia_maxReIm  find the maximum of all real and of all imaginary parts separately along the diagonal of a square matrix; merge these two maxima into one complex number and return it as a scalar 
MCF_Rows_sabsmax  find the complex numbers of largest sum Re+Im along rows and store these row maxima in a (complex) column vector 
MCF_Cols_sabsmax  find the complex numbers of largest sum Re+Im along columns and store these column maxima in a (complex) row vector 
MCF_Dia_sabsmax  find the complex number of largest sum Re+Im along the diagonal of a square matrix and return it as a (complex) scalar 
To determine the center of gravity of a matrix, you have the choice between the following two functions:
MF_centerOfGravityInd  center of gravity, returned as interpolated element indices 
MF_centerOfGravityV  center of gravity of an MZ matrix with explicitly given X and Y axes 
Back to MatrixLib Table of Contents OptiVec home
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 
Back to MatrixLib Table of Contents OptiVec home
MF_addM  add two matrices 
MF_addMT  add one matrix and the transpose of another matrix MC = MA + MB^{T} 
MF_subM  subtract one matrix from another 
MF_subMT  subtract a transposed matrix MC = MA  MB^{T} 
MF_subrMT  subtract a matrix from another, transposed, matrix MC = MB^{T}  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 + MB^{T}) 
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  MB^{T}) 
MFs_subrMT  scaled reverse subtraction of one matrix and the transpose of another: MC = c * (MB^{T}  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 = MA^{T} * 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 * MA^{T} 
MF_mulM  multiply two matrices: MC = MA * MB 
MF_mulMT  multiply one matrix by the transpose of another matrix: MC = MA * MB^{T} 
MF_TmulM  multiply the transpose of a matrix by another matrix: MC = MA^{T} * MB 
MF_TmulMT  multiply the transpose of one matrix by the transpose of another matrix: MC = MA^{T} * MB^{T} 
MCF_mulMH  multiply one matrix by the hermitian conjugate of another matrix: MC = MA * MB^{T *} 
MCF_HmulM  multiply the hermitian conjugate of a matrix by another matrix: MC = MA^{T *} * 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 = MA^{T} * 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 * MB^{T} 
Back to MatrixLib Table of Contents OptiVec home
MF_solve  solve simultaneous linear equations (using LU decomposition) 
MF_inv  invert a matrix 
MF_det  determinant of a matrix 
MF_solveBySVD  solve simultaneous linear equations, using Singular Value Decomposition 
MF_safeSolve  tries first solution by LUD; if that fails, SVD is done 
MF_LUdecompose  decompose into LU form 
MF_LUimprove  improve accuracy of LU decomposition by iteration 
MF_LUDresult  check if MF_LUdecompose was successful 
MF_LUDsetEdit  set editing threshold for MF_LUdecompose; may be used to work around singularities 
MF_LUDgetEdit  retrieve currently set threshold 
MF_LUsolve  solve simultaneous linear equations, given the matrix in LU form 
MF_LUinv  invert matrix already composed into LU form 
MF_LUdet  determinant of matrix already composed into LU form 
MF_SVdecompose  Singular Value Decomposition 
MF_SVsolve  solve SVdecomposed set of linear equations 
MF_SVDsetEdit  set threshold for Singular Value editing 
MF_SVDgetEdit  retrieve current threshold for Singular Value editing 
Back to MatrixLib Table of Contents OptiVec home
MFsym_eigenvalues  eigenvalues with or without eigenvectors of a symmetric real matrix 
MFsym_sqrt  squareroot of a symmetric, positive definite matrix 
Back to MatrixLib Table of Contents OptiVec home
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 crosscorrelation 
MF_spectrum  Spatial frequency spectrum 
MF_Rows_FFT  Fourier Transform along rows 
MF_Cols_FFT  Fourier Transform along columns 
Back to MatrixLib Table of Contents OptiVec home
VF_linregress  equallyweighted linear regression on XY data 
VF_linregresswW  the same with nonequal weighting 
VF_polyfit  fitting of one XY data set to a polynomial 
VF_polyfitwW  the same for nonequal datapoint weighting 
VF_linfit  fitting of one XY data set to an arbitrary function linear in its parameters 
VF_linfitwW  the same for nonequal datapoint weighting 
VF_nonlinfit  fitting of one XY data set to an arbitrary, possibly nonlinear function 
VF_nonlinfitwW  the same for nonequal datapoint weighting 
VF_multiLinfit  fitting of multiple XY data sets to one common linear function 
VF_multiLinfitwW  the same for nonequal datapoint weighting 
VF_multiNonlinfit  fitting of multiple XY data sets to one common nonlinear function 
VF_multiNonlinfitwW  the same for nonequal datapoint weighting 
MF_linfit  fit XYZ data to a function linear in its parameters 
MF_linfitwW  the same with nonequal weighting of individual data points 
MF_nonlinfit  fit XYZ data to an arbitrary, possibly nonlinear function 
MF_nonlinfitwW  the same for nonequal datapoint weighting 
MF_multiLinfit  fit multiple XYZ data sets to one common linear function 
MF_multiLinfitwW  the same for nonequal datapoint weighting 
MF_multiNonlinfit  fit multiple XYZ data sets to one common nonlinear function 
MF_multiNonlinfitwW  the same for nonequal datapoint weighting 
void LinModel( fVector BasFuncs, float x, unsigned nfuncs )
{
BasFuncs[0] = sin(x);
BasFuncs[1] = cos(x);
}
Note that the coefficients a_{i} are not used in the model function itself. The argument nfuncs (which is neglected in the above example) allows to use a variable number of basis functions. It is also possible to switch fitting parameters "on" and "off", i.e., "free" or "fixed". To this end, the routine VF_linfit takes a "status" array as an argument. For all members of the parameter vector that are to be treated as free, the corresponding "status" entry must be set to 1. If status_{i} is set to 0, the corresponding parameter, a_{i}, is treated as fixed at its value upon input.
Internally, VF_linfit employs a Singular Value Decomposition algorithm to obtain a solution even for (near)singular linear systems. Thereby, coefficients a_{i} whose significance is lower than a threshold, Thresh, are
set to 0 instead of infinity. This threshold can be modified by calling VF_setLinfitNeglect. The current threshold can be retrieved by VF_getLinfitNeglect.
Having noted above that functions like
y = sin(a_{0}x)+ cos(a_{1}x)
require their own treatment, we arrive at the last and most general class of model functions:
void NonlinModel( fVector Y, fVector X, ui size)
{
for( ui i=0; i<size; i++ )
Y[i] = sin( a[0]*X[i] ) + cos( a[1]*X[i] );
}
Keeping in mind that a nonlinear fitting routine will spend most of its time in your model function (and its derivatives, see below), you might wish to optimize this function, 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)
{
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" should be global. It is passed as an argument to VF_nonlinfit, and its elements are changed during the fitting process. 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.
All nonlinfit functions return the figureofmerit of the best parameter array "a" found. To perform the fit, these functions need not only the fitting function itself, but also its partial derivatives with respect to the individual coefficients. Therefore, an argument "Derivatives" is required in calling the nonlinfit functions. If you happen to know the partial derivatives analytically, Derivatives should point to a function calculating them. If you do not know them, call with Derivatives = NULL (nil in Pascal/Delphi). In the latter case, all derivatives will be calculated numerically inside the functions. In cases where you know some, but not all of the partial derivatives of Y with respect to the coefficients a_{i}, you could also calculate dY / da_{i} wherever you have an analytic formula, and call VF_nonlinfit_autoDeriv for the remaining coefficients. To demonstrate this possibility, here is a function coding the derivatives for our nonlinear sine and cosine example:
VF_NONLINFITWORKSPACE ws; // &ws must be passed both to VF_nonlinfit and to VF_nonlinfit_autoDeriv
void DerivModel( fVector dYdAi, fVector X, ui size, unsigned iPar )
{
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);
}
}
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 5! 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 )
{
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 nonlinear 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
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.
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 c^{2} (default: EPSILON) float FracTolChi; fractional change of c^{2} (default: SQRT_EPSILON) float AbsTolPar; absolute change of all parameters (default: SQRT_MIN) float FracTolPar; fractional change of all parameters (default: SQRT_EPSILON)
unsigned HowOftenFulfill; how often fulfill one of the above conditions before convergence is considered achieved (default: 3) unsigned LevelOfMethod; 1: LevenbergMarquardt 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 c^{2} (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 reshaping 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 userdefined 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.
As described above, all nonlinear 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 userretrievable 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.
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) 
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) 
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.
VF_nonlinfit_getBestValues  best set of parameters found so far 
VF_nonlinfit_getChi2  the best figureofmerit (c^{2}, chisquare) obtained so far 
VF_nonlinfit_getTestDir  returns the test direction (+1 for upwards, 1 for downwards) during "breakout" attempts (levelofmethod 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 LevenbergMarquardt or Downhill cycle 
Back to MatrixLib Table of Contents OptiVec home
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 cutoff); (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 
Back to MatrixLib Table of Contents OptiVec home
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 colordensity plots. 
M_setDensityMapBounds  Set a color scale and draw an XY coordinate system for matrix colordensity plots. 
M_findDensityMapBounds  Calculate a color scale and draw an XY coordinate system for matrix colordensity 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 © 19982013 OptiCode  Dr. Martin Sander Software Development