# NAG Library Routine Document

## 1Purpose

e04rnf is a part of the NAG optimization modelling suite and defines one or more linear matrix constraints of the problem.

## 2Specification

Fortran Interface
 Subroutine e04rnf ( nvar, dima, nnza, a, nblk,
 Integer, Intent (In) :: nvar, dima, nnza(nvar+1), nnzasum, irowa(nnzasum), icola(nnzasum), nblk, blksizea(nblk) Integer, Intent (Inout) :: idblk, ifail Real (Kind=nag_wp), Intent (In) :: a(nnzasum) Type (c_ptr), Intent (In) :: handle
#include nagmk26.h
 void e04rnf_ (void **handle, const Integer *nvar, const Integer *dima, const Integer nnza[], const Integer *nnzasum, const Integer irowa[], const Integer icola[], const double a[], const Integer *nblk, const Integer blksizea[], Integer *idblk, Integer *ifail)

## 3Description

After the initialization routine e04raf has been called, e04rnf may be used to add one or more linear matrix inequalities
 $∑ i=1 n xi Ai - A0 ⪰ 0$ (1)
to the problem definition. Here ${A}_{i}$ are $d$ by $d$ symmetric matrices. The expression $S⪰0$ stands for a constraint on eigenvalues of a symmetric matrix $S$, namely, all the eigenvalues should be non-negative, i.e., the matrix $S$ should be positive semidefinite.
Typically, this will be used in linear semidefinite programming problems (SDP)
 $minimize x∈ℝn cTx (a) subject to ∑ i=1 n xi Aik - A0k ⪰ 0 , k=1,…,mA (b) lB≤Bx≤uB (c) lx≤x≤ux (d)$ (2)
or to define the linear part of bilinear matrix inequalities (3)(b) in (BMI-SDP)
 $minimize x∈ℝn 12 xTHx + cTx (a) subject to ∑ i,j=1 n xi xj Qijk + ∑ i=1 n xi Aik - A0k ⪰ 0 , k=1,…,mA (b) lB≤Bx≤uB (c) lx≤x≤ux (d)$ (3)
e04rnf can be called repeatedly to accumulate more matrix inequalities. See e04raf for more details.

### 3.1Input data organization

All the matrices ${A}_{\mathit{i}}$, for $\mathit{i}=0,1,\dots ,n$, are symmetric and thus only their upper triangles are passed to the routine. They are stored in sparse coordinate storage format (see Section 2.1.1 in the F11 Chapter Introduction), i.e., every nonzero from the upper triangles is coded as a triplet of row index, column index and the numeric value. These triplets of all (upper triangle) nonzeros from all ${A}_{i}$ matrices are passed to the routine in three arrays: irowa for row indices, icola for column indices and a for the values. No particular order of nonzeros within one matrix is enforced but all nonzeros from ${A}_{0}$ must be stored first, followed by all nonzero from ${A}_{1}$, followed by ${A}_{2}$, etc.
The number of stored nonzeros from each ${A}_{i}$ matrix is given in ${\mathbf{nnza}}\left(i+1\right)$, thus this array indicates which section of arrays irowa, icola and a belongs to which ${A}_{i}$ matrix. See Table 1 and the example in Section 9. See also e04rdf which uses the same data organization.
 irowa upper triangle upper triangle upper triangle icola nonzeros nonzeros $\cdots$ nonzeros a $\underbrace{\text{from ​}{A}_{0}}$ $\underbrace{\text{from ​}{A}_{1}}$ $\underbrace{\text{from ​}{A}_{n}}$ ${\mathbf{nnza}}\left(1\right)$ ${\mathbf{nnza}}\left(2\right)$ ${\mathbf{nnza}}\left(n+1\right)$
Table 1
Coordinate storage format of matrices ${A}_{0},{A}_{1},\dots ,{A}_{n}$ in input arrays
There are two possibilities for defining more matrix inequality constraints
 $∑ i=1 n xi A i k - A 0 k ⪰ 0 , k=1,2,…,mA$ (4)
to the problem. The first is to call e04rnf ${m}_{A}$ times and define a single matrix inequality at a time. This might be more straightforward and therefore it is recommended. Alternatively, it is possible to merge all ${m}_{A}$ constraints into one inequality and pass them in a single call to e04rnf. It is easy to see that (4) can be equivalently expressed as one bigger matrix inequality with the following block diagonal structure
 $∑ i=1 n xi Ai1 Ai2 ⋱ AimA - A01 A02 ⋱ A0mA ⪰ 0 .$
If ${d}_{k}$ denotes the dimension of inequality $k$, the new merged inequality has dimension $d=\sum _{\mathit{k}=1}^{{m}_{A}}{d}_{k}$ and each of the ${A}_{i}$ matrices is formed by ${A}_{i}^{1},{A}_{i}^{2},\dots ,{A}_{i}^{{m}_{A}}$ stored as ${m}_{A}$ diagonal blocks. In such a case, nblk is set to ${m}_{A}$ and ${\mathbf{blksizea}}\left(k\right)$ to ${d}_{k}$, the size of the $k$th diagonal blocks. This might be useful in connection with e04rdf.
On the other hand, if there is no block structure and just one matrix inequality is provided, nblk should be set to $1$ and blksizea is not referenced.

### 3.2Definition of Bilinear Matrix Inequalities (BMI)

e04rnf is designed to be used together with e04rpf to define bilinear matrix inequalities (3)(b). e04rnf sets the linear part of the constraint and e04rpf expands it by higher order terms. To distinquish which linear matrix inequality (or more precisely, which block) is to be expanded, e04rpf needs the number of the block, idblk. The blocks are numbered as they are added, starting from $1$.
Whenever a matrix inequality (or a set of them expressed as diagonal blocks) is stored, the routine returns idblk of the last inequality added. idblk is just the order of the inequality amongst all matrix inequalities accumulated through the calls. The first inequality has ${\mathbf{idblk}}=1$, the second one ${\mathbf{idblk}}=2$, etc. Therefore if you call e04rnf for the very first time with ${\mathbf{nblk}}=42$, it adds $42$ inequalities with idblk from $1$ to $42$ and the routine returns ${\mathbf{idblk}}=42$ (the number of the last one). A subsequent call with ${\mathbf{nblk}}=1$ would add only one inequality, this time with ${\mathbf{idblk}}=43$ which would be returned.

None.

## 5Arguments

1:     $\mathbf{handle}$ – Type (c_ptr)Input
On entry: the handle to the problem. It needs to be initialized by e04raf and must not be changed.
2:     $\mathbf{nvar}$ – IntegerInput
On entry: $n$, the number of decision variables $x$ in the problem. It must be unchanged from the value set during the initialization of the handle by e04raf.
3:     $\mathbf{dima}$ – IntegerInput
On entry: $d$, the dimension of the matrices ${A}_{\mathit{i}}$, for $\mathit{i}=0,1,\dots ,{\mathbf{nvar}}$.
Constraint: ${\mathbf{dima}}>0$.
4:     $\mathbf{nnza}\left({\mathbf{nvar}}+1\right)$ – Integer arrayInput
On entry: ${\mathbf{nnza}}\left(\mathit{i}+1\right)$, for $\mathit{i}=0,1,\dots ,{\mathbf{nvar}}$, gives the number of nonzero elements in the upper triangle of matrix ${A}_{i}$. To define ${A}_{i}$ as a zero matrix, set ${\mathbf{nnza}}\left(i+1\right)=0$. However, there must be at least one matrix with at least one nonzero.
Constraints:
• ${\mathbf{nnza}}\left(i\right)\ge 0$;
• $\sum _{\mathit{i}=1}^{n+1}{\mathbf{nnza}}\left(i\right)\ge 1$.
5:     $\mathbf{nnzasum}$ – IntegerInput
On entry: the dimension of the arrays irowa, icola and a, at least the total number of all nonzeros in all matrices ${A}_{i}$.
Constraints:
• ${\mathbf{nnzasum}}>0$;
• $\sum _{\mathit{i}=1}^{n+1}{\mathbf{nnza}}\left(i\right)\le {\mathbf{nnzasum}}$.
6:     $\mathbf{irowa}\left({\mathbf{nnzasum}}\right)$ – Integer arrayInput
7:     $\mathbf{icola}\left({\mathbf{nnzasum}}\right)$ – Integer arrayInput
8:     $\mathbf{a}\left({\mathbf{nnzasum}}\right)$ – Real (Kind=nag_wp) arrayInput
On entry: nonzero elements in upper triangle of matrices ${A}_{i}$ stored in coordinate storage. The first ${\mathbf{nnza}}\left(1\right)$ elements belong to ${A}_{0}$, the following ${\mathbf{nnza}}\left(2\right)$ elements belong to ${A}_{1}$, etc. See explanation above.
Constraints:
• $1\le {\mathbf{irowa}}\left(i\right)\le {\mathbf{dima}}$, ${\mathbf{irowa}}\left(i\right)\le {\mathbf{icola}}\left(i\right)\le {\mathbf{dima}}$;
• irowa and icola match the block diagonal pattern set by blksizea.
9:     $\mathbf{nblk}$ – IntegerInput
On entry: ${m}_{A}$, number of diagonal blocks in ${A}_{i}$ matrices. As explained above it is equivalent to the number of matrix inequalities supplied in this call.
Constraint: ${\mathbf{nblk}}\ge 1$.
10:   $\mathbf{blksizea}\left({\mathbf{nblk}}\right)$ – Integer arrayInput
On entry: if ${\mathbf{nblk}}>1$, sizes ${d}_{k}$ of the diagonal blocks.
If ${\mathbf{nblk}}=1$, blksizea is not referenced.
Constraints:
• ${\mathbf{blksizea}}\left(i\right)\ge 1$;
• $\sum _{\mathit{i}=1}^{{m}_{A}}{\mathbf{blksizea}}\left(i\right)={\mathbf{dima}}$.
11:   $\mathbf{idblk}$ – IntegerInput/Output
On entry: if ${\mathbf{idblk}}=0$, new matrix inequalities are created. This is the only value allowed at the moment; nonzero values are reserved for future releases of the NAG Library.
Constraint: ${\mathbf{idblk}}=0$.
On exit: the number of the last matrix inequality added. By definition, it is the number of the matrix inequalities already defined plus nblk.
12:   $\mathbf{ifail}$ – IntegerInput/Output
On entry: ifail must be set to $0$, $-1\text{​ or ​}1$. If you are unfamiliar with this argument you should refer to Section 3.4 in How to Use the NAG Library and its Documentation for details.
For environments where it might be inappropriate to halt program execution when an error is detected, the value $-1\text{​ or ​}1$ is recommended. If the output of error messages is undesirable, then the value $1$ is recommended. Otherwise, the recommended value is $-1$. When the value $-\mathbf{1}\text{​ or ​}1$ is used it is essential to test the value of ifail on exit.
On exit: ${\mathbf{ifail}}={\mathbf{0}}$ unless the routine detects an error or a warning has been flagged (see Section 6).

## 6Error Indicators and Warnings

If on entry ${\mathbf{ifail}}=0$ or $-1$, explanatory error messages are output on the current error message unit (as defined by x04aaf).
Errors or warnings detected by the routine:
${\mathbf{ifail}}=1$
The supplied handle does not define a valid handle to the data structure for the NAG optimization modelling suite. It has not been initialized by e04raf or it has been corrupted.
${\mathbf{ifail}}=2$
The problem cannot be modified in this phase any more, the solver has already been called.
${\mathbf{ifail}}=4$
On entry, ${\mathbf{idblk}}=〈\mathit{\text{value}}〉$.
Constraint: ${\mathbf{idblk}}=0$.
On entry, ${\mathbf{nvar}}=〈\mathit{\text{value}}〉$, expected $\mathrm{value}=〈\mathit{\text{value}}〉$.
Constraint: nvar must match the value given during initialization of handle.
${\mathbf{ifail}}=6$
On entry, ${\mathbf{dima}}=〈\mathit{\text{value}}〉$.
Constraint: ${\mathbf{dima}}>0$.
On entry, $i=〈\mathit{\text{value}}〉$ and ${\mathbf{nnza}}\left(i\right)=〈\mathit{\text{value}}〉$.
Constraint: ${\mathbf{nnza}}\left(i\right)\ge 0$.
On entry, ${\mathbf{nnzasum}}=〈\mathit{\text{value}}〉$ and $\mathrm{sum}\left({\mathbf{nnza}}\right)=〈\mathit{\text{value}}〉$.
Constraint: ${\mathbf{nnzasum}}\ge \mathrm{sum}\left({\mathbf{nnza}}\right)$.
On entry, $\mathrm{sum}\left({\mathbf{nnza}}\right)=〈\mathit{\text{value}}〉$.
Constraint: $\mathrm{sum}\left({\mathbf{nnza}}\right)\ge 1$.
${\mathbf{ifail}}=7$
On entry, ${\mathbf{dima}}=〈\mathit{\text{value}}〉$ and $\mathrm{sum}\left({\mathbf{blksizea}}\right)=〈\mathit{\text{value}}〉$.
Constraint: $\mathrm{sum}\left({\mathbf{blksizea}}\right)={\mathbf{dima}}$.
On entry, $i=〈\mathit{\text{value}}〉$ and ${\mathbf{blksizea}}\left(i\right)=〈\mathit{\text{value}}〉$.
Constraint: ${\mathbf{blksizea}}\left(i\right)\ge 1$.
On entry, ${\mathbf{nblk}}=〈\mathit{\text{value}}〉$.
Constraint: ${\mathbf{nblk}}>0$.
${\mathbf{ifail}}=8$
An error occurred in matrix ${A}_{i}$, $i=〈\mathit{\text{value}}〉$ (counting indices $1\dots {\mathbf{nvar}}+1$).
On entry, $j=〈\mathit{\text{value}}〉$, ${\mathbf{icola}}\left(j\right)=〈\mathit{\text{value}}〉$ and ${\mathbf{dima}}=〈\mathit{\text{value}}〉$.
Constraint: $1\le {\mathbf{icola}}\left(j\right)\le {\mathbf{dima}}$.
An error occurred in matrix ${A}_{i}$, $i=〈\mathit{\text{value}}〉$ (counting indices $1\dots {\mathbf{nvar}}+1$).
On entry, $j=〈\mathit{\text{value}}〉$, ${\mathbf{irowa}}\left(j\right)=〈\mathit{\text{value}}〉$ and ${\mathbf{dima}}=〈\mathit{\text{value}}〉$.
Constraint: $1\le {\mathbf{irowa}}\left(j\right)\le {\mathbf{dima}}$.
An error occurred in matrix ${A}_{i}$, $i=〈\mathit{\text{value}}〉$ (counting indices $1\dots {\mathbf{nvar}}+1$).
On entry, $j=〈\mathit{\text{value}}〉$, ${\mathbf{irowa}}\left(j\right)=〈\mathit{\text{value}}〉$ and ${\mathbf{icola}}\left(j\right)=〈\mathit{\text{value}}〉$.
Constraint: ${\mathbf{irowa}}\left(\mathit{j}\right)\le {\mathbf{icola}}\left(\mathit{j}\right)$ (elements within the upper triangle).
An error occurred in matrix ${A}_{i}$, $i=〈\mathit{\text{value}}〉$ (counting indices $1\dots {\mathbf{nvar}}+1$).
On entry, $j=〈\mathit{\text{value}}〉$, ${\mathbf{irowa}}\left(j\right)=〈\mathit{\text{value}}〉$ and ${\mathbf{icola}}\left(j\right)=〈\mathit{\text{value}}〉$. Maximum column index in this row given by the block structure defined by blksizea is $〈\mathit{\text{value}}〉$.
Constraint: all elements of ${A}_{i}$ must respect the block structure given by blksizea.
An error occurred in matrix ${A}_{i}$, $i=〈\mathit{\text{value}}〉$ (counting indices $1\dots {\mathbf{nvar}}+1$).
On entry, more than one element of ${A}_{i}$ has row index $〈\mathit{\text{value}}〉$ and column index $〈\mathit{\text{value}}〉$.
Constraint: each element of ${A}_{i}$ must have a unique row and column index.
${\mathbf{ifail}}=-99$
See Section 3.9 in How to Use the NAG Library and its Documentation for further information.
${\mathbf{ifail}}=-399$
Your licence key may have expired or may not have been installed correctly.
See Section 3.8 in How to Use the NAG Library and its Documentation for further information.
${\mathbf{ifail}}=-999$
Dynamic memory allocation failed.
See Section 3.7 in How to Use the NAG Library and its Documentation for further information.

Not applicable.

## 8Parallelism and Performance

e04rnf is not threaded in any implementation.

The following example demonstrates how the elements of the ${A}_{i}^{k}$ matrices are organized within the input arrays. Let us assume that there are two blocks defined (${\mathbf{nblk}}=2$). The first has dimension $3$ by $3$ (${\mathbf{blksizea}}\left(1\right)=3$) and the second $2$ by $2$ (${\mathbf{blksizea}}\left(2\right)=2$). For simplicity, the number of variables is $2$. Please note that the values were chosen to ease orientation rather than to define a valid problem.
 $A 0 1 = 0.1 0 0.3 0 0.2 0.4 0.3 0.4 0 , A 1 1 ​ empty ​ A 2 1 = 2.1 0 0 0 2.2 0 0 0 2.3 , A 0 2 = 0 -0.1 -0.1 0 , A 1 2 = -1.1 0 0 -1.2 , A 2 2 = -2.1 -2.2 -2.2 -2.3 .$
Both inequalities will be passed in a single call to e04rnf, therefore the matrices are merged into the following block diagonal form:
 $A0 = 0.1 0 0.3 0 0.2 0.4 0.3 0.4 0 0 -0.1 -0.1 0 ,$
 $A1 = 0 0 0 0 0 0 0 0 0 -1.1 0 0 -1.2 ,$
 $A2 = 2.1 0 0 0 2.2 0 0 0 2.3 -2.1 -2.2 -2.2 -2.3 .$
All matrices are symmetric and therefore only the upper triangles are passed to the routine. The coordinate storage format is used. Note that elements within the same matrix do not need to be in any specific order. The table below shows one of the ways the arrays could be populated.
 irowa $\begin{array}{ccccc}\phantom{0.}2\hfill & \phantom{0.}2\hfill & \phantom{-0.}4\hfill & \phantom{0.}1\hfill & \phantom{0.}1\hfill \end{array}$ $\begin{array}{cc}\phantom{-0.}4& \phantom{-0.}5\end{array}$ $\begin{array}{cccccc}\phantom{0.}1& \phantom{0.}2& \phantom{0.}3& \phantom{-0.}4& \phantom{-0.}4& \phantom{-0.}5\end{array}$ icola $\begin{array}{ccccc}\phantom{0.}2& \phantom{0.}3& \phantom{-0.}5& \phantom{0.}1& \phantom{0.}3\end{array}$ $\begin{array}{cc}\phantom{-0.}4& \phantom{-0.}5\end{array}$ $\begin{array}{cccccc}\phantom{0.}1& \phantom{0.}2& \phantom{0.}3& \phantom{-0.}4& \phantom{-0.}5& \phantom{-0.}5\end{array}$ a $\underbrace{\begin{array}{ccccc}0.2& 0.4& -0.1& 0.1& 0.3\end{array}}$ $\underbrace{\begin{array}{cc}-1.1& -1.2\end{array}}$ $\underbrace{\begin{array}{cccccc}2.1& 2.2& 2.3& -2.1& -2.2& -2.3\end{array}}$ ${A}_{0}$ ${A}_{1}$ ${A}_{2}$ nnza $5$ $2$ $6$

## 10Example

There are various problems which can be successfully reformulated and solved as an SDP problem. The following example shows how a maximization of the minimal eigenvalue of a matrix depending on certain parameters can be utilized in statistics.
For further examples, please refer to Section 10 in e04raf.
Given a series of $M$ vectors of length $p$, $\left\{{v}_{i}:i=1,2,\dots ,M\right\}$ this example solves the SDP problem:
 $maximize λ1,…,λM,t t subject to ∑ i=1 M λi vi viT ⪰ t⁢I ∑ i=1 M λi =1 λi≥0 , k=1,…,M .$
This formulation comes from an area of statistics called experimental design and corresponds to finding an approximate $E$ optimal design for a linear regression.
A linear regression model has the form:
 $y=Xβ+ε$
where $y$ is a vector of observed values, $X$ is a design matrix of (known) independent variables and $\epsilon$ is a vector of errors. In experimental design it is assumed that each row of $X$ is chosen from a set of $M$ possible vectors, $\left\{{v}_{i}:i=1,2,\dots ,M\right\}$. The goal of experimental design is to choose the rows of $X$ so that the error covariance is ‘small’. For an $E$ optimal design this is defined as the $X$ that maximizes the minimum eigenvalue of ${X}^{\mathrm{T}}X$.
In this example we construct the $E$ optimal design for a polynomial regression model of the form:
 $y = β0+ β1⁢ x+ β2⁢ x2+ β3⁢ x3+ β4⁢ x4+ε$
where $x\in \left\{1-j×0.05:j=0,1,\dots ,40\right\}$.

### 10.1Program Text

Program Text (e04rnfe.f90)

### 10.2Program Data

Program Data (e04rnfe.d)

### 10.3Program Results

Program Results (e04rnfe.r)

© The Numerical Algorithms Group Ltd, Oxford, UK. 2017