Solution to a Generalized Linear Regression Model Problem

using a generalized QR factorization of A and B, where A is an N×M matrix, B is an N×P matrix, and d is a vector of length N. It is also assumed that M N M+P and:

Under these assumptions, the constrained equation is always consistent, and there is a unique solution x and a minimal 2-norm solution y. In particular, if matrix B is square nonsingular, then the problem GLM is equivalent to the following weighted linear least squares problem:

Calling Sequence

CALL DGGGLM	(N, M, P, DA, LDA, DB, LDB, DD, DX, DY, DWORK, LDWORK, INFO)
CALL SGGGLM	(N, M, P, SA, LDA, SB, LDB, SD, SX, SY, SWORK, LDWORK, INFO)
CALL ZGGGLM	(N, M, P, ZA, LDA, ZB, LDB, ZD, ZX, ZY, ZWORK, LDWORK, INFO)
CALL CGGGLM	(N, M, P, CA, LDA, CB, LDB, CD, CX, CY, CWORK, LDWORK, INFO)
void dggglm	(int n, int m, int p, double *da, int lda, double *db, int ldb, double *dd, double *dx, double *dy, int *info)
void sggglm	(int n, int m, int p, float *sa, int lda, float *sb, int ldb, float *sd, float *sx, float *dy, int *info)
void zggglm	(int n, int m, int p, doublecomplex *za, int lda, doublecomplex *zb, int ldb, doublecomplex *zd, doublecomplex *zx, doublecomplex *zy, int *info)
void cggglm	(int n, int m, int p, complex *ca, int lda, complex *cb, int ldb, complex *cd, complex *cx, complex *cy, int *info)

Arguments

N	Number of rows of the matrices A and B. N 0.
M	Number of columns of the matrix A. M 0.
P	Number of columns of the matrix B. M+P N M.
xA	On entry, the N×M matrix A. On exit, A is overwritten.
LDA	Leading dimension of the array A as specified in a dimension or type statement. LDA max(1, N).
xB	On entry, the N×P matrix B. On exit, B is overwritten.
LDB	Leading dimension of the array B as specified in a dimension or type statement. LDB max(1, N).
xD	On entry, the left-hand side of the GLM equation. On exit, D is overwritten.
xX	On exit, X and Y are the solutions of the GLM problem.
xY	On exit, X and Y are the solutions of the GLM problem.
xWORK	Scratch array with a dimension of LDWORK.
LDWORK	Leading dimension of the array WORK as specified in a dimension or type statement. LDWORK M + P + max(N, M, P).
INFO	On exit:
	INFO = 0	Subroutine completed normally.
	INFO < 0	The ith argument, where i = |INFO|, had an illegal value.

Sample Program

PROGRAM TEST

IMPLICIT NONE

C

INTEGER LDA, LDB, LDWORK, N

PARAMETER (N = 3)

PARAMETER (LDA = N)

PARAMETER (LDB = N)

PARAMETER (LDWORK = 3 * N)

C

DOUBLE PRECISION A(LDA,N), B(LDB,N), DLEFT(N), X(N), Y(N)

DOUBLE PRECISION WORK(LDWORK)

INTEGER ICOL, INFO, IROW

C

EXTERNAL DGGGLM

INTRINSIC ABS

C

C Initialize the array A to store the matrix A shown below.

C Initialize the array B to store the matrix B shown below.

C Initialize the array DLEFT to store the left hand side d

C of the GLM equation as shown below.

C

C 1 4 6 11 -8 2 1

C A = 2 9 14 B = -4 5 -2 d = 10

C 3 14 23 1 -2 1 100

C

DATA A / 1.0D0, 2.0D0, 3.0D0,

$ 4.0D0, 9.0D0, 1.4D1,

$ 6.0D0, 1.4D1, 2.3D1 /

DATA B / 1.1D1, -4.0D0, 1.0D0,

$ -8.0D0, 5.0D0, 2.0D0,

$ 2.0D0, -2.0D0, 1.0D0 /

DATA DLEFT / 1.0D0, 1.0D1, 1.0D2 /

C

C Print the initial values of A and B.

C

PRINT 1000

PRINT 1010, ((A(IROW,ICOL), ICOL = 1, N), IROW = 1, N)

PRINT 1020

PRINT 1010, ((B(IROW,ICOL), ICOL = 1, N), IROW = 1, N)

C

C Compute and print the solution to the linear regression

C problem.

C

CALL DGGGLM (N, N, N, A, LDA, B, LDB, DLEFT, X, Y, WORK,

$ LDWORK, INFO)

IF (INFO .LT. 0) THEN

PRINT 1030, ABS(INFO)

STOP 1

END IF

PRINT 1040

DO 100, IROW = 1, N

PRINT 1050, X(IROW), Y(IROW)

100 CONTINUE

C

1000 FORMAT (1X, 'A:')

1010 FORMAT (3(3X, F4.1))

1020 FORMAT (/1X, 'B:')

1030 FORMAT (1X, 'Illegal value for argument #', I1,

$ ' in DGGGLM.')

1040 FORMAT (/1X, 3X, 'X', 7X, 'Y')

1050 FORMAT (1X, F6.1, 2X, F6.1)

C

END

Sample Output

A:
1.0 4.0 6.0
2.0 9.0 14.0
3.0 14.0 23.0
B:
11.0 -8.0 2.0
-4.0 5.0 -2.0
1.0 2.0 1.0
X Y
131.0 0.0
-154.0 0.0
81.0 0.0