Matrix functions take one or more matrices (or attributes) as arguments. Unlike most other LINGO functions, which primarily return scalar (or single-value) results, some matrix functions can return multiple results, with some of the results being matrices or vectors.  Matrix functions supported are:

 

L[, err] = @CHOLESKY( A)

The @CHOLESKY function performs a Cholesky factorization of the symmetric, positive-definite matrix A, where LL' = A, and L is a lower-triangular matrix with real and positive diagonal elements.  One can think of L as the "square root" of matrix A.  Both A and L must be square and of the same dimension.  The upper triangle of L will be filled with zeros.

The optional argument, err, is an error flag that will be 0 if the factorization was successful, otherwise it will be 1.  If err is present, any error will be ignored by LINGO and will be left for your model to handle. If err is not present, then, on any error, LINGO will display an error message and terminate the run.

In this example below, we take a small 3x3 matrix, A, compute its Cholesky factor, L, then show that LL' is equal to the original matrix A.

MODEL:

! Compute the Cholesky factorization of matrix A.

! Back check by taking the Cholesky factor, L, and

! multiplying it with its transpose, L', to get

! the original matrix A;

SETS:
 S1;

 S2(S1,S1): A, A2, L, LT;

ENDSETS

 

DATA:

  S1 = 1..3;

  A = 10 11 12

      11 54 51

      12 51 86;

ENDDATA

 

CALC:

  ! Get the Cholesky factor;

  L = @CHOLESKY( A);

 

  ! L transpose;

  LT = @TRANSPOSE( L);

 

  ! Compute L * L';

  @FOR( S2( I, J):

      A2( I, J) = @SUM( S1( J2): L( I, J2) * LT( J2, J));

  );

ENDCALC

 

CALC:

  @SET( 'TERSEO', 2)  ; ! Just display errors;

  @WRITE( "Original A matrix:", @NEWLINE( 1));

  @TABLE( A);

  @WRITE( @NEWLINE( 1), "Cholesky factor L:", @NEWLINE( 1));

  @TABLE( L);

  @WRITE( @NEWLINE( 1), "L transpose:", @NEWLINE( 1));

  @TABLE( LT);

  @WRITE( @NEWLINE( 1), "L * L' = A:", @NEWLINE( 1));

  @TABLE( A2);

ENDCALC

END

Model: CHOLESKY

The output from this model is reproduced below, where we can see, as per definition, LL' = A:

Original A matrix:

            1         2         3

  1  10.00000  11.00000  12.00000

  2  11.00000  54.00000  51.00000

  3  12.00000  51.00000  86.00000

 

Cholesky factor L:

            1         2         3

  1  3.162278  0.000000  0.000000

  2  3.478505  6.473021  0.000000

  3  3.794733  5.839623  6.123627

 

L transpose:

            1         2         3

  1  3.162278  3.478505  3.794733

  2  0.000000  6.473021  5.839623

  3  0.000000  0.000000  6.123627

 

L * L' = A:

            1         2         3

  1  10.00000  11.00000  12.00000

  2  11.00000  54.00000  51.00000

  3  12.00000  51.00000  86.00000

 

A = @MTXMUL( B, C)

The @MTXMUL function multiplies matrix B by matrix C and places the result in matrix A.  The matrices must all be defined on dense sets.  The dimensions of the matrices must also agree, for example, if B is an m x n matrix, then C must be an n x p, and A must be an m x p matrix.  Furthermore, m, n and p must be greater than 1, i.e., inner and outer products are not currently supported.

In the Cholesky example above, you will recall that we computed L x L', i.e., the product of matrix L by its transpose, L'.  We explicitly wrote out the formula for matrix multiplication with the following code:

  ! Compute L * L';

  @FOR( S2( I, J):

      A2( I, J) = @SUM( S1( J2): L( I, J2) * LT( J2, J));

  );

A simpler and more efficient way to have done this would have been to use the @MTXMUL function as follows:

  ! Compute L * L';

  A2 = @MTXMUL( L, LT);

 

@DETERMINANT( A)

The @DETERMINANT function returns the determinant value for matrix A.

In the example below, we make use of Cramer's Rule for finding the solution to a system of linear equations of the form Ax=b.  According to Cramer's Rule, the solution for x( i) may be calculated by computing the determinant of the A matrix with column i replaced by the right-hand side values b. This determinant value is then divided by the determinant of the original A matrix to get the final value for x( i).

MODEL:

 

! Use Cramer's Rule to find the solution

 to a system of linear equations, Ax=b;

 

SETS:

  DIM: X, B, B2;

  DXD( DIM, DIM): A, A2;

ENDSETS

 

DATA:

  DIM = 1..3;

  B = 12 0 -3;

  A =

     1  -2   5

     4   1  -2

    -2   1  -1

  ;

ENDDATA

 

CALC:

  ! Determinant of the original matrix;

  D = @DETERMINANT( A);

 

  ! Determine if a unique solution exists;

  @IFC( D #EQ# 0:

     @STOP( '  No unique solution exists!!!');

  @ELSE

     @WRITE( '  A unique solution exists...');

  );

  @WRITE( @NEWLINE( 1));

 

  ! Copy matrix;

  @FOR( DXD( I, J): A2( I, J) = A( I, J));

 

  @FOR( DIM( J):

     ! Replace column J with the B values;        

     @FOR( DIM( I):

        A2( I, J) = B( I);

     );

 

     ! Calculate solution for X( J);

     X( J) = @DETERMINANT( A2) / D;

 

     ! Restore column J in the original matrix;

     @FOR( DIM( I):

        A2( I, J) = A( I, J);

     );

  );

 

  ! Back check the solution;  

  @FOR( DIM( I):

     B2( I) = @SUM( DIM( J): A( I, J) * X( J));

     @IFC( @ABS( B2( I) - B( I)) #GT# 1.E-6:

        @STOP( 'Solution fails!!!!');

     );

  );

  @WRITE( '  Solution verified...', @NEWLINE( 1),

   '  Solution = ');

 

  @FOR( DIM( I): @WRITE( X( I), '   '));

 

  @WRITE( @NEWLINE( 2));

ENDCALC

 

END

Model: CRAMER

We first checked to see if a unique solution to the system exists by computing the determinant of the original matrix and checking that it was not 0.  We then used Cramer's Rule to compute the solution vector, x. Finally, we took the x values and verified that they solve the equations.  The output from the model is:

A unique solution exists...

Solution verified...

Solution = 1   2   3

 
LAMDAR, VR, LAMBDAI, VI, ERR = @EIGEN( A)

The @EIGEN function returns eigenvectors and eigenvalues for the matrix A, where:

As long as at least one left-hand side argument is present, all other arguments may be omitted.  If err is present, then LINGO will not halt the run if a numeric error occurs; it will be up to your model to handle the error.  If err is not present and a numeric error occurs, then LINGO will halt model execution.

In the next example, we compute the eigenvalues and eigenvectors for a 3x3 symmetric matrix.

MODEL:

 

SETS:

  DIM: LAMBDA;

  DXD( DIM, DIM): A, V, AV, LV;

  DXDL( DXD) | &1 #GT# &2: PROD;

ENDSETS

 

DATA:

  TOL = 1.E-6;

  DIM = 1..3;

  A =  3 4 1

       4 0 9

       1 9 5;

ENDDATA

 

CALC:

  ! Compute the eigenvalues and eigen vectors;

  ! Note: The matrix is symmetric, so the vectors

  ! and values will not have imaginary parts, so

  ! we can omit the imaginary arguments.;

  LAMBDA, V,,, ERR = @EIGEN( A);

 

  ! Numeric error?;

  @IFC( ERR #NE# 0: @STOP( 'Numeric error'));

 
  ! Display them with the @TABLE function.;

  @WRITE( @NEWLINE( 1), ' Eigenvalues:', @NEWLINE( 2));

  @TABLE( LAMBDA);

 

  @WRITE( @NEWLINE( 1), ' Eigenvectors:', @NEWLINE( 1));

  @TABLE( V);

 

  ! Check that A * v = Lambda * v holds.;

  NZ = 0;

  @FOR( DXD( I, J):

     LV( I, J) = LAMBDA( J) * V( I, J);

     AV( I, J) = @SUM( DIM( K): A( I, K) * V( K, J));

     @IFC( @ABS( LV( I, J) - AV( I, J)) #GT# TOL: NZ = NZ + 1);

  );

  @WRITE( @NEWLINE( 1));

  @IFC( NZ #GT# 0:

     @WRITE( ' Lambda * v = A * v is not verified.');

  @ELSE

     @WRITE( ' Lambda * v = A * v is verified.');

  );

 
  ! Symmetric matrices have orthogonal eigenvectors.

  ! Verify this by computing the inner products.;

  NZ = 0;

  @FOR( DXDL( I, J):

     PROD( I, J) = @sum( DIM( K): V( K, I) * V( K, J));

     @IFC( @ABS( PROD( I, J)) #GT# TOL: NZ = NZ + 1);

  );

  @WRITE( @NEWLINE( 1));

  @IFC( NZ #GT# 0:

     @WRITE( ' Eigenvectors are NOT orthogonal.', @NEWLINE( 2));

  @ELSE

     @WRITE( ' Eigenvectors are orthogonal.', @NEWLINE( 2));

  );

ENDCALC

END

Model: EIGENEX1

Symmetric matrices have the property that their eigenvalues and eigenvectors are all real valued, allowing us to ignore the imaginary components.  Looking at the call to @EIGEN:

LAMBDA, V,,, ERR = @EIGEN( A);

we see that two left-hand side arguments are omitted, corresponding to the imaginary parts.  Note, we could have passed these two arguments, however, they would have simply been filled with all zeros.

By the definition of eigenvalues and eigenvectors, it's true that Av = Lambda * v, where A is the original nxn matrix, v is an nx1 eigenvector of A, and Lambda is the eigenvalue corresponding to v. The model verifies this relationship and displays the message:

Lambda * v = A * v is verified.

Another nice property of symmetric matrices is that the eigenvectors are all orthogonal.  The model verifies this property by calculating all the inner products between the eigenvectors and comparing these products to 0.

The output from the model follows:

Eigenvalues:

 

 1   12.92039

 2   2.587500

 3  -7.507889

 

Eigenvectors:

                 1               2               3

 1      -0.3179176      -0.9145159       0.2501781

 2      -0.6062189  -0.6813624E-02      -0.7952686

 3      -0.7289904       0.4044926       0.5522307

 

Lambda * v = A * v is verified.

Eigenvectors are orthogonal.

 

AINV[, err] = @INVERSE( A)

The @INVERSE function returns the inverse of matrix A.  Both AINV and A must be square, nxn matrices.

The optional argument, err, is an error flag that will be 0 if the invert operation was successful, otherwise it will be 1.  A situation where err would be nonzero is when A is a singular matrix.  If err is present, any error will be ignored by LINGO and will be left for your model to handle. If err is not present, then, on any error, LINGO will display an error message and terminate the run.

The example below may be familiar from the @DETERMINANT example above, where we solve a 3x3 system of linear equations. In this example, however, we will not use Cramer's Rule to solve the system.  Instead, we will compute the inverse of the A matrix, and then multiply the inverse by the right-hand side vector b to get the solution vector x. We also verify that a unique solution exists by computing the determinant of A, as well as verify that our solution does indeed solve the original set of equations.

MODEL:

 

SETS:

  DIM: X, B, B2;

  DXD( DIM, DIM): A, AINV;

ENDSETS

 

DATA:

  DIM = 1..3;

  B = 12 0 -3;

  A =

     1  -2   5

     4   1  -2

    -2   1  -1

  ;

ENDDATA

 

CALC:

  ! Determinant of the original matrix;

  D = @DETERMINANT( A);

 

  ! Determine if a unique solution exists;

  @IFC( D #EQ# 0:

     @STOP( '  No unique solution exists!!!');

  @ELSE

     @WRITE( '  A unique solution exists...');

  );

  @WRITE( @NEWLINE( 1));

 

  ! Get the matrix inverse;

  AINV = @INVERSE( A);

 

  ! Use the inverse to compute the solution vector: x = AINV * b;

  @FOR( DIM( I): X( I) = @SUM( DIM( J): AINV( I, J) * B( J)));

 

  ! Back check the solution;  

  @FOR( DIM( I):

     B2( I) = @SUM( DIM( J): A( I, J) * X( J));

     @IFC( @ABS( B2( I) - B( I)) #GT# 1.E-6:

        @STOP( 'Solution fails!!!!');

     );

  );

  @WRITE( '  Solution verified...', @NEWLINE( 1),

   '  Solution = ');

 

  @FOR( DIM( I): @WRITE( @FORMAT( X( I), '10.5g')));

 

  @WRITE( @NEWLINE( 2));

ENDCALC

 

END

Model: INVERSE

The output from the model follows:

A unique solution exists...

Solution verified...

Solution =          1         2         3

Again, as with the @DETERMINANT example, the solution is [1,2,3].

 

B, b0, RES, rsq, f, p, var = @REGRESS( Y, X)

The @REGRESS function is used to perform multiple linear regression (MLR), a technique that models the linear relationship, Y = b0 + B*X, between one dependent variable and one or more independent variables.  MLR is frequently used in forecasting, e.g., by economists in macroeconomics when modeling of the economy.

Using the syntax above, we have the following definitions:

Note, that any of the left-hand side arguments may be omitted if they are not of interest.  The only requirement is that at least one left-hand side argument must be present.

The goal of MLR is to fit a linear function to the observed data, X, for predicting the independent observations, Y. The form of the function is:

 Y = b0 + ∑ i B( i) * X( i) + RES

The fit is determined by the set of coefficients that minimizes the sum of the squared residual values.  For more information, see: https://en.wikipedia.org/wiki/Line.

In the LINGO sample models folder, there is a regression example called REGRLONG, shown below. This model makes use of the Longley dataset:

http://www.itl.nist.gov/div898/strd/lls/data/LINKS/DATA/Longley.dat

This is a series of seven macroeconomic statistics, with one being total employment, the dependent variable.  The Longley dataset is commonly used for testing statistical software, and the link above contains the certified results for regressions based on the data.

MODEL:

 

! Multiple linear regression;

! The output is the set of regression coefficients, B,

for the model:

   Y(i) = B(0) + B(1)*X(i,1) + B(2)*X(i,2)+... + ERROR(i)

 

We use the Longley macroeconomic data set for:

http://www.itl.nist.gov/div898/strd/lls/data/LINKS/DATA/Longley.dat

 

Certified Regression Statistics for Longley data:

                                        Standard Deviation

 Parameter         Estimate                      of Estimate

  B(0)          -3482258.63459582          890420.383607373

  B(1)           15.0618722713733          84.9149257747669

  B(2)          -0.358191792925910E-01     0.334910077722432E-01

  B(3)          -2.02022980381683          0.488399681651699

  B(4)          -1.03322686717359          0.214274163161675

  B(5)          -0.511041056535807E-01     0.226073200069370

  B(6)           1829.15146461355          455.478499142212

 

Residual

Standard Deviation   304.854073561965                

 

R-Squared            0.995479004577296                

 

Certified Analysis of Variance Table:

 

Source of   Degrees of       Sums of              Mean

Variation    Freedom          Squares             Squares

Regression      6      184172401.944494    30695400.3240823

Residual        9      836424.055505915    92936.0061673238

 

  F Statistic

330.285339234588

;

 

SETS:

 OBS : Y, RES;

 VARS;

 DEPVAR( VARS);    ! The dependent variable;

 EXPVAR( VARS): B; ! The explanatory variables;

 OXV( OBS, VARS): DATATAB; ! The data table;

 OXE( OBS, EXPVAR): X;

ENDSETS

 

DATA:

 

! Names of the variables;

 VARS = EMPLD PRICE GNP__ JOBLS MILIT POPLN  YEAR_;

 

! Dependent variable. Must be exactly 1;

 DEPVAR = EMPLD;

 

! Explanatory variables. Should not include

  dependent variable;

 EXPVAR = PRICE GNP__ JOBLS MILIT POPLN YEAR_;

 

! The dataset;

 DATATAB =

   60323   83  234289  2356  1590  107608  1947

   61122  88.5 259426  2325  1456  108632  1948

   60171  88.2 258054  3682  1616  109773  1949

   61187  89.5 284599  3351  1650  110929  1950

   63221  96.2 328975  2099  3099  112075  1951

   63639  98.1 346999  1932  3594  113270  1952

   64989  99   365385  1870  3547  115094  1953

   63761 100   363112  3578  3350  116219  1954

   66019 101.2 397469  2904  3048  117388  1955

   67857 104.6 419180  2822  2857  118734  1956

   68169 108.4 442769  2936  2798  120445  1957

   66513 110.8 444546  4681  2637  121950  1958

   68655 112.6 482704  3813  2552  123366  1959

   69564 114.2 502601  3931  2514  125368  1960

   69331 115.7 518173  4806  2572  127852  1961

   70551 116.9 554894  4007  2827  130081  1962;

 

ENDDATA

 

CALC:

   !Suppress some of Lingo's output;

   @SET( 'TERSEO', 2);

 

   !Set up data for the @REGRESS function

    by moving the dependent variable into

    Y and the independents into X;

   @FOR( DEPVAR( K):

      @FOR( OBS( I):

         Y( I) = DATATAB( I, K);

         @FOR( OXE( I, J): X( I, J) = DATATAB( I, J))

      );

   );

 

   ! Do the regression;

   B, B0, RES, RSQUARE, F, P, VAR = @REGRESS( Y, X);

 

   ! Print a report;

   @WRITE( '   Explained error/R-Square: ',

    @FORMAT( RSQUARE, '16.8g'), @NEWLINE(1));

   @WRITE( '                F statistic: ',

    @FORMAT( F, '16.8g'), @NEWLINE(1));

   @WRITE( '                    p value: ',

    @FORMAT( P, '16.8g'), @NEWLINE(1));

   @WRITE( '          Residual variance: ',

    @FORMAT( VAR, '16.8g'), @NEWLINE(1));

   @WRITE( '         Residual std error: ',

    @FORMAT( VAR ^ .5, '16.8g'), @NEWLINE(2));

   @WRITE( '   B( 0):     ', @FORMAT( B0, '16.8g'),

    @NEWLINE( 1));

   @FOR( EXPVAR( I):

    @WRITE( '   B( ', EXPVAR( I), '): ',

    @FORMAT( B( I), '16.8g'), @NEWLINE( 1));

   );

ENDCALC

END

Model: REGRLONG

The interesting parts of this model are the following two expressions:

   !Set up data for the @REGRESS function

    by moving the dependent variable into

    Y and the independents into X;

   @FOR( OBS( I):

      Y( I) = DATATAB( I, @INDEX( EMPLD));

      @FOR( OXE( I, J): X( I, J) = DATATAB( I, J))

   );

 

   ! Do the regression;

   B, B0, RES, RSQUARE, F, P, VAR = @REGRESS( Y, X);

Note, that the data is provided as a single table, containing both the dependent and independent observations.  The first expression partitions the raw data into a vector Y of dependent employment observations, and a 16x6 X matrix of independent observations.  At this point, the data is in the correct format to submit to @REGRESS, and we perform the regression with the statement:

    ! Do the regression;

B, B0, RES, RSQUARE, F, P, VAR = @REGRESS( Y, X);

Results from the regression are then displayed:

  Explained error/R-Square:         0.995479

               F statistic:        330.28534

                   p value:    4.984031e-010

         Residual variance:        92936.006

        Residual std error:        304.85407

 

  B( 0):           -3482258.6

  B( PRICE):        15.061872

  B( GNP__):     -0.035819179

  B( JOBLS):       -2.0202298

  B( MILIT):       -1.0332269

  B( POPLN):     -0.051104106

  B( YEAR_):        1829.1515

 

T = @TRANSPOSE( A)

The @TRANSPOSE function returns the transpose of matrix A, where T(i,j) = A(j,i) for every (i,j).  Matrix A does not have to be square, but if A is an m x n matrix, then T must be of dimension nxm.  For an example of @TRANSPOSE, refer to the section above on @CHOLESKY, where we use @TRANSPOSE to get the transpose of the Cholesky factor.  This transpose is then multiplied by the original Cholesky factor to get the original matrix back.

Q,R[, err] = @QRFACTOR( A)

The @QRFACTOR function returns a QR factorization of matrix A, where A = QR, with R being upper triangular and Q being orthogonal (i.e., Q'Q is the identity matrix).  The matrices must all be defined on dense sets, with A and Q being m x n matrices and R an n x n matrix.

The optional argument, err, is an error flag that will be 0 if the factorization operation was successful, otherwise it will be 1.  If err is present, any error will be ignored by LINGO and will be left for your model to handle. If err is not present, then, on any error, LINGO will display an error message and terminate the run.

The sample model QRFACTOR contained in the LINGO Samples folder illustrates the use of the @QRFACTOR function.