# 9.1 Solving Linear Systems Involving the Basis Matrix¶

A linear optimization problem always has an optimal solution which is also a basic solution. In an optimal basic solution there are exactly \(m\) basic variables where \(m\) is the number of rows in the constraint matrix \(A\). Define

as a matrix consisting of the columns of \(A\) corresponding to the basic variables. The basis matrix \(B\) is always non-singular, i.e.

or, equivalently, \(B^{-1}\) exists. This implies that the linear systems

and

each have a unique solution for all \(w\).

**MOSEK** provides functions for solving the linear systems (9.1) and (9.2) for an arbitrary \(w\).

In the next sections we will show how to use **MOSEK** to

## 9.1.1 Basis identification¶

To use the solutions to (9.1) and (9.2) it is important to know how the basis matrix \(B\) is constructed.

Internally **MOSEK** employs the linear optimization problem

where

The basis matrix is constructed of \(m\) columns taken from

If variable \(x_j\) is a basis variable, then the \(j\)-th column of \(A\), denoted \(a_{:,j}\), will appear in \(B\). Similarly, if \(x_i^c\) is a basis variable, then the \(i\)-th column of \(-I\) will appear in the basis. The ordering of the basis variables and therefore the ordering of the columns of \(B\) is arbitrary. The ordering of the basis variables may be retrieved by calling the function `initbasissolve`

. This function initializes data structures for later use and returns the indexes of the basic variables in the array `basis`

. The interpretation of the `basis`

is as follows. If we have

then the \(i\)-th basis variable is

Moreover, the \(i\)-th column in \(B\) will be the \(i\)-th column of \(-I\). On the other hand if

then the \(i\)-th basis variable is the variable

and the \(i\)-th column of \(B\) is the column

For instance if \(\mathtt{basis}[0] = 4\) and \(\mathtt{numcon} = 5\), then since \(\mathtt{basis}[0]< \mathtt{numcon}\), the first basis variable is \(x_4^c\). Therefore, the first column of \(B\) is the fourth column of \(-I\). Similarly, if \(\mathtt{basis}[1] = 7\), then the second variable in the basis is \(x_{\mathtt{basis}[1]-\mathtt{numcon}} = x_2\). Hence, the second column of \(B\) is identical to \(a_{:,2}\).

An example

Consider the linear optimization problem:

Suppose a call to `initbasissolve`

returns an array `basis`

so that

```
basis[0] = 1,
basis[1] = 2.
```

Then the basis variables are \(x_1^c\) and \(x_0\) and the corresponding basis matrix \(B\) is

Please note the ordering of the columns in \(B\) .

```
using Mosek
maketask() do task
# Use remote server: putoptserverhost(task,"http://solve.mosek.com:30080")
putobjname(task,"solvebasis")
putstreamfunc(task,MSK_STREAM_LOG,msg -> print(msg))
numcon = 2
numvar = 2
c = Float64[1.0, 1.0]
ptrb = Int64[1,3]
ptre = Int64[3,4]
asub = Int32[1,2,
1,2]
aval = Float64[1.0, 1.0,
2.0, 1.0]
bkc = Boundkey[MSK_BK_UP
MSK_BK_UP]
blc = Float64[-Inf
-Inf]
buc = Float64[2.0
6.0]
bkx = Boundkey[MSK_BK_LO
MSK_BK_LO]
blx = Float64[0.0
0.0]
bux = Float64[Inf
Inf]
w1 = Float64[2.0, 6.0]
w2 = Float64[1.0, 0.0]
inputdata(task,
Int32(numcon), Int32(numvar),
c,
0.0,
ptrb,
ptre,
asub,
aval,
bkc,
blc,
buc,
bkx,
blx,
bux)
putobjsense(task,MSK_OBJECTIVE_SENSE_MAXIMIZE)
r = optimize(task)
if r != MSK_RES_OK
println("Mosek warning: $r")
end
basis = initbasissolve(task)
#List basis variables corresponding to columns of B
varsub = Int32[1, 2]
for i in 1:numcon
if basis[varsub[i]] < numcon
println("Basis variable no $i is xc$(basis[i])")
else
println("Basis variable no $i is x$(basis[i]-numcon)")
# solve Bx = w1
# varsub contains index of non-zeros in b.
# On return b contains the solution x and
# varsub the index of the non-zeros in x.
nz = 2
nz = solvewithbasis(task,false, nz, varsub, w1)
println("nz = $nz")
println("Solution to Bx = $w1")
for i in 1:nz
if basis[varsub[i]] < numcon
println("xc $(basis[varsub[i]]) = $(w1[varsub[i]])")
else
println("x$(basis[varsub[i]] - numcon) = $(w1[varsub[i]])")
end
end
# Solve B^Tx = w2
nz = 1
varsub[1] = 1
nz = solvewithbasis(task,true,nz,varsub,w2)
println(nz)
println("Solution to B^Tx = $(w2)")
for i in 1:nz
if basis[varsub[i]] < numcon
println("xc$(basis[varsub[i]]) = $(w2[varsub[i]])")
else
println("x$(basis[varsub[i]] - numcon) = $(w2[varsub[i]])")
end
end
end
end
end
```

In the example above the linear system is solved using the optimal basis for (9.4) and the original right-hand side of the problem. Thus the solution to the linear system is the optimal solution to the problem. When running the example program the following output is produced.

```
basis[0] = 1
Basis variable no 0 is xc1.
basis[1] = 2
Basis variable no 1 is x0.
Solution to Bx = b:
x0 = 2.000000e+00
xc1 = -4.000000e+00
Solution to B^Tx = c:
x1 = -1.000000e+00
x0 = 1.000000e+00
```

Please note that the ordering of the basis variables is

and thus the basis is given by:

It can be verified that

is a solution to

## 9.1.2 Solving arbitrary linear systems¶

**MOSEK** can be used to solve an arbitrary (rectangular) linear system

using the `solvewithbasis`

function without optimizing the problem as in the previous example. This is done by setting up an \(A\) matrix in the task, setting all variables to basic and calling the `solvewithbasis`

function with the \(b\) vector as input. The solution is returned by the function.

An example

Below we demonstrate how to solve the linear system

with two inputs \(b=(1,-2)\) and \(b=(7,0)\) .

```
using Mosek
function setup(task :: Mosek.Task,
aval :: Vector{Float64},
asub :: Vector{Int32},
ptrb :: Vector{Int64},
ptre :: Vector{Int64},
numvar :: Int32)
appendvars(task,numvar)
appendcons(task,numvar)
putacolslice(task,1,numvar+1,ptrb,ptre,asub,aval)
putconboundsliceconst(task,1,numvar+1,MSK_BK_FX,0.0,0.0)
putvarboundsliceconst(task,1,numvar+1,MSK_BK_FR,-Inf,Inf)
# Define a basic solution by specifying status keys for variables
# & constraints.
deletesolution(task,MSK_SOL_BAS)
putskcslice(task,MSK_SOL_BAS, 1, numvar+1, fill(MSK_SK_FIX,numvar))
putskxslice(task,MSK_SOL_BAS, 1, numvar+1, fill(MSK_SK_BAS,numvar))
return initbasissolve(task)
end
let numcon = Int32(2),
numvar = Int32(2),
aval = [ -1.0 ,
1.0, 1.0 ],
asub = Int32[ 2,
1, 2 ],
ptrb = Int64[1, 2],
ptre = Int64[2, 4]
# bsub = new int[numvar];
# b = new double[numvar];
# basis = new int[numvar];
maketask() do task
# Use remote server: putoptserverhost(task,"http://solve.mosek.com:30080")
putstreamfunc(task,MSK_STREAM_LOG,msg -> print(msg))
# Put A matrix and factor A. Call this function only once for a
# given task.
basis = setup(task,
aval,
asub,
ptrb,
ptre,
numvar)
# now solve rhs
let b = Float64[1,-2],
bsub = Int32[1,2],
nz = solvewithbasis(task,false, 2, bsub, b)
println("Solution to Bx = b:")
# Print solution and show correspondents to original variables in the problem
for i in 1:nz
if basis[bsub[i]] <= numcon
println("This should never happen")
else
println("x $(basis[bsub[i]] - numcon) = $(b[bsub[i]])")
end
end
end
let b = Float64[7,0],
bsub = Int32[1,0],
nz = solvewithbasis(task,false,1, bsub, b)
println("Solution to Bx = b:")
# Print solution and show correspondents to original variables in the problem
for i in 1:nz
if (basis[bsub[i]] <= numcon)
println("This should never happen")
else
println("x $(basis[bsub[i]] - numcon) = $(b[bsub[i]])")
end
end
end
end
end
```

The most important step in the above example is the definition of the basic solution, where we define the status key for each variable. The actual values of the variables are not important and can be selected arbitrarily, so we set them to zero. All variables corresponding to columns in the linear system we want to solve are set to basic and the slack variables for the constraints, which are all non-basic, are set to their bound.

The program produces the output:

```
Solution to Bx = b:
x1 = 1
x0 = 3
Solution to Bx = b:
x1 = 7
x0 = 7
```