text:
# Interactive Simplex Method¶

This module, meant for educational purposes only, supports learning and exploring of the simplex method.

Do you want to solve Linear Programs efficiently? use MixedIntegerLinearProgram instead.

The methods implemented here allow solving Linear Programming Problems (LPPs) in a number of ways, may require explicit (and correct!) description of steps and are likely to be much slower than “regular” LP solvers. If, however, you want to learn how the simplex method works and see what happens in different situations using different strategies, but don’t want to deal with tedious arithmetic, this module is for you!

Historically it was created to complement the Math 373 course on Mathematical Programming and Optimization at the University of Alberta, Edmonton, Canada.

AUTHORS:

• Andrey Novoseltsev (2013-03-16): initial version.
• Matthias Koeppe, Peijun Xiao (2015-07-05): allow different output styles.

EXAMPLES:

Most of the module functionality is demonstrated on the following problem.

Corn & Barley

A farmer has 1000 acres available to grow corn and barley. Corn has a net profit of 10 dollars per acre while barley has a net profit of 5 dollars per acre. The farmer has 1500 kg of fertilizer available with 3 kg per acre needed for corn and 1 kg per acre needed for barley. The farmer wants to maximize profit. (Sometimes we also add one more constraint to make the initial dictionary infeasible: the farmer has to use at least 40% of the available land.)

Using variables $$C$$ and $$B$$ for land used to grow corn and barley respectively, in acres, we can construct the following LP problem:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblem(A, b, c, ["C", "B"], variable_type=">=")
sage: P
LP problem (use typeset mode to see details)

It is recommended to copy-paste such examples into your own worksheet, so that you can run these commands with typeset mode on and get

$\begin{split}\begin{array}{l} \begin{array}{lcrcrcl} \max \mspace{-6mu}&\mspace{-6mu} \mspace{-6mu}&\mspace{-6mu} 10 C \mspace{-6mu}&\mspace{-6mu} + \mspace{-6mu}&\mspace{-6mu} 5 B \mspace{-6mu} \\ \mspace{-6mu}&\mspace{-6mu} \mspace{-6mu}&\mspace{-6mu} C \mspace{-6mu}&\mspace{-6mu} + \mspace{-6mu}&\mspace{-6mu} B \mspace{-6mu}&\mspace{-6mu} \leq \mspace{-6mu}&\mspace{-6mu} 1000 \\ \mspace{-6mu}&\mspace{-6mu} \mspace{-6mu}&\mspace{-6mu} 3 C \mspace{-6mu}&\mspace{-6mu} + \mspace{-6mu}&\mspace{-6mu} B \mspace{-6mu}&\mspace{-6mu} \leq \mspace{-6mu}&\mspace{-6mu} 1500 \\ \end{array} \\ C, B \geq 0 \end{array}\end{split}$

Since it has only two variables, we can solve it graphically:

sage: P.plot()
Graphics object consisting of 19 graphics primitives

The simplex method can be applied only to problems in standard form, which can be created either directly

sage: InteractiveLPProblemStandardForm(A, b, c, ["C", "B"])
LP problem (use typeset mode to see details)

or from an already constructed problem of “general type”:

sage: P = P.standard_form()

In this case the problem does not require any modifications to be written in standard form, but this step is still necessary to enable methods related to the simplex method.

The simplest way to use the simplex method is:

sage: P.run_simplex_method()
\begin{equation*}
...
The optimal value: $6250$. An optimal solution: $\left(250,\,750\right)$.

(This method produces quite long formulas which have been omitted here.) But, of course, it is much more fun to do most of the steps by hand. Let’s start by creating the initial dictionary:

sage: D = P.initial_dictionary()
sage: D
LP problem dictionary (use typeset mode to see details)

Using typeset mode as recommended, you’ll see

$\begin{split}\renewcommand{\arraystretch}{1.5} \begin{array}{|rcrcrcr|} \hline x_{3} \mspace{-6mu}&\mspace{-6mu} = \mspace{-6mu}&\mspace{-6mu} 1000 \mspace{-6mu}&\mspace{-6mu} - \mspace{-6mu}&\mspace{-6mu} C \mspace{-6mu}&\mspace{-6mu} - \mspace{-6mu}&\mspace{-6mu} B\\ x_{4} \mspace{-6mu}&\mspace{-6mu} = \mspace{-6mu}&\mspace{-6mu} 1500 \mspace{-6mu}&\mspace{-6mu} - \mspace{-6mu}&\mspace{-6mu} 3 C \mspace{-6mu}&\mspace{-6mu} - \mspace{-6mu}&\mspace{-6mu} B\\ \hline z \mspace{-6mu}&\mspace{-6mu} = \mspace{-6mu}&\mspace{-6mu} 0 \mspace{-6mu}&\mspace{-6mu} + \mspace{-6mu}&\mspace{-6mu} 10 C \mspace{-6mu}&\mspace{-6mu} + \mspace{-6mu}&\mspace{-6mu} 5 B\\ \hline \end{array}\end{split}$

With the initial or any other dictionary you can perform a number of checks:

sage: D.is_feasible()
True
sage: D.is_optimal()
False

You can look at many of its pieces and associated data:

sage: D.basic_variables()
(x3, x4)
sage: D.basic_solution()
(0, 0)
sage: D.objective_value()
0

Most importantly, you can perform steps of the simplex method by picking an entering variable, a leaving variable, and updating the dictionary:

sage: D.enter("C")
sage: D.leave(4)
sage: D.update()

If everything was done correctly, the new dictionary is still feasible and the objective value did not decrease:

sage: D.is_feasible()
True
sage: D.objective_value()
5000

If you are unsure about picking entering and leaving variables, you can use helper methods that will try their best to tell you what are your next options:

sage: D.possible_entering()
[B]
sage: D.possible_leaving()
Traceback (most recent call last):
...
ValueError: leaving variables can be determined
for feasible dictionaries with a set entering variable
or for dual feasible dictionaries

It is also possible to obtain feasible sets and final dictionaries of problems, work with revised dictionaries, and use the dual simplex method!

Note

Currently this does not have a display format for the terminal.

## Classes and functions¶

class sage.numerical.interactive_simplex_method.InteractiveLPProblem(A, b, c, x='x', constraint_type='<=', variable_type='', problem_type='max', base_ring=None, is_primal=True, objective_constant_term=0)

Construct an LP (Linear Programming) problem.

Note

This class is for educational purposes only: if you want to solve Linear Programs efficiently, use MixedIntegerLinearProgram instead.

This class supports LP problems with “variables on the left” constraints.

INPUT:

• A – a matrix of constraint coefficients
• b – a vector of constraint constant terms
• c – a vector of objective coefficients
• x – (default: "x") a vector of decision variables or a string giving the base name
• constraint_type – (default: "<=") a string specifying constraint type(s): either "<=", ">=", "==", or a list of them
• variable_type – (default: "") a string specifying variable type(s): either ">=", "<=", "" (the empty string), or a list of them, corresponding, respectively, to non-negative, non-positive, and free variables
• problem_type – (default: "max") a string specifying the problem type: "max", "min", "-max", or "-min"
• base_ring – (default: the fraction field of a common ring for all input coefficients) a field to which all input coefficients will be converted
• is_primal – (default: True) whether this problem is primal or dual: each problem is of course dual to its own dual, this flag is mostly for internal use and affects default variable names only
• objective_constant_term – (default: 0) a constant term of the objective

EXAMPLES:

We will construct the following problem:

$\begin{split}\begin{array}{l} \begin{array}{lcrcrcl} \max \mspace{-6mu}&\mspace{-6mu} \mspace{-6mu}&\mspace{-6mu} 10 C \mspace{-6mu}&\mspace{-6mu} + \mspace{-6mu}&\mspace{-6mu} 5 B \mspace{-6mu} \\ \mspace{-6mu}&\mspace{-6mu} \mspace{-6mu}&\mspace{-6mu} C \mspace{-6mu}&\mspace{-6mu} + \mspace{-6mu}&\mspace{-6mu} B \mspace{-6mu}&\mspace{-6mu} \leq \mspace{-6mu}&\mspace{-6mu} 1000 \\ \mspace{-6mu}&\mspace{-6mu} \mspace{-6mu}&\mspace{-6mu} 3 C \mspace{-6mu}&\mspace{-6mu} + \mspace{-6mu}&\mspace{-6mu} B \mspace{-6mu}&\mspace{-6mu} \leq \mspace{-6mu}&\mspace{-6mu} 1500 \\ \end{array} \\ C, B \geq 0 \end{array}\end{split}$
sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblem(A, b, c, ["C", "B"], variable_type=">=")

Same problem, but more explicitly:

sage: P = InteractiveLPProblem(A, b, c, ["C", "B"],
....:     constraint_type="<=", variable_type=">=")

Even more explicitly:

sage: P = InteractiveLPProblem(A, b, c, ["C", "B"], problem_type="max",
....:     constraint_type=["<=", "<="], variable_type=[">=", ">="])

Using the last form you should be able to represent any LP problem, as long as all like terms are collected and in constraints variables and constants are on different sides.

A()

Return coefficients of constraints of self, i.e. $$A$$.

OUTPUT:

• a matrix

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblem(A, b, c, ["C", "B"], variable_type=">=")
sage: P.constraint_coefficients()
[1 1]
[3 1]
sage: P.A()
[1 1]
[3 1]

Abcx()

Return $$A$$, $$b$$, $$c$$, and $$x$$ of self as a tuple.

OUTPUT:

• a tuple

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblem(A, b, c, ["C", "B"], variable_type=">=")
sage: P.Abcx()
(
[1 1]
[3 1], (1000, 1500), (10, 5), (C, B)
)

add_constraint(coefficients, constant_term, constraint_type='<=')

Return a new LP problem by adding a constraint toself.

INPUT:

• coefficients – coefficients of the new constraint
• constant_term – a constant term of the new constraint
• constraint_type – (default: "<=") a string indicating the constraint type of the new constraint

OUTPUT:

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblem(A, b, c)
sage: P1 = P.add_constraint(([2, 4]), 2000, "<=")
sage: P1.Abcx()
(
[1 1]
[3 1]
[2 4], (1000, 1500, 2000), (10, 5), (x1, x2)
)
sage: P1.constraint_types()
('<=', '<=', '<=')
sage: P.Abcx()
(
[1 1]
[3 1], (1000, 1500), (10, 5), (x1, x2)
)
sage: P.constraint_types()
('<=', '<=')
sage: P2 = P.add_constraint(([2, 4, 6]), 2000, "<=")
Traceback (most recent call last):
...
TypeError: number of columns must be the same, not 2 and 3
sage: P3 = P.add_constraint(([2, 4]), 2000, "<")
Traceback (most recent call last):
...
ValueError: unknown constraint type

b()

Return constant terms of constraints of self, i.e. $$b$$.

OUTPUT:

• a vector

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblem(A, b, c, ["C", "B"], variable_type=">=")
sage: P.constant_terms()
(1000, 1500)
sage: P.b()
(1000, 1500)

base_ring()

Return the base ring of self.

Note

The base ring of LP problems is always a field.

OUTPUT:

• a ring

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblem(A, b, c, ["C", "B"], variable_type=">=")
sage: P.base_ring()
Rational Field

sage: c = (10, 5.)
sage: P = InteractiveLPProblem(A, b, c, ["C", "B"], variable_type=">=")
sage: P.base_ring()
Real Field with 53 bits of precision

c()

Return coefficients of the objective of self, i.e. $$c$$.

OUTPUT:

• a vector

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblem(A, b, c, ["C", "B"], variable_type=">=")
sage: P.objective_coefficients()
(10, 5)
sage: P.c()
(10, 5)

constant_terms()

Return constant terms of constraints of self, i.e. $$b$$.

OUTPUT:

• a vector

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblem(A, b, c, ["C", "B"], variable_type=">=")
sage: P.constant_terms()
(1000, 1500)
sage: P.b()
(1000, 1500)

constraint_coefficients()

Return coefficients of constraints of self, i.e. $$A$$.

OUTPUT:

• a matrix

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblem(A, b, c, ["C", "B"], variable_type=">=")
sage: P.constraint_coefficients()
[1 1]
[3 1]
sage: P.A()
[1 1]
[3 1]

constraint_types()

Return a tuple listing the constraint types of all rows.

OUTPUT:

• a tuple of strings

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblem(A, b, c, ["C", "B"], variable_type=">=", constraint_type=["<=", "=="])
sage: P.constraint_types()
('<=', '==')

decision_variables()

Return decision variables of self, i.e. $$x$$.

OUTPUT:

• a vector

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblem(A, b, c, ["C", "B"], variable_type=">=")
sage: P.decision_variables()
(C, B)
sage: P.x()
(C, B)

dual(y=None)

Construct the dual LP problem for self.

INPUT:

• y – (default: depends on style()) a vector of dual decision variables or a string giving the base name

OUTPUT:

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblem(A, b, c, ["C", "B"], variable_type=">=")
sage: DP = P.dual()
sage: DP.b() == P.c()
True
sage: DP.dual(["C", "B"]) == P
True

feasible_set()

Return the feasible set of self.

OUTPUT:

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblem(A, b, c, ["C", "B"], variable_type=">=")
sage: P.feasible_set()
A 2-dimensional polyhedron in QQ^2
defined as the convex hull of 4 vertices

is_bounded()

Check if self is bounded.

OUTPUT:

• True is self is bounded, False otherwise

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblem(A, b, c, ["C", "B"], variable_type=">=")
sage: P.is_bounded()
True

Note that infeasible problems are always bounded:

sage: b = (-1000, 1500)
sage: P = InteractiveLPProblem(A, b, c, variable_type=">=")
sage: P.is_feasible()
False
sage: P.is_bounded()
True

is_feasible(*x)

Check if self or given solution is feasible.

INPUT:

• (optional) anything that can be interpreted as a valid solution for this problem, i.e. a sequence of values for all decision variables

OUTPUT:

• True is this problem or given solution is feasible, False otherwise

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblem(A, b, c, variable_type=">=")
sage: P.is_feasible()
True
sage: P.is_feasible(100, 200)
True
sage: P.is_feasible(1000, 200)
False
sage: P.is_feasible([1000, 200])
False
sage: P.is_feasible(1000)
Traceback (most recent call last):
...
TypeError: given input is not a solution for this problem

is_negative()

Return $$True$$ when the problem is of type "-max" or "-min".

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblem(A, b, c, ["C", "B"], variable_type=">=")
sage: P.is_negative()
False
sage: P = InteractiveLPProblem(A, b, c, ["C", "B"], variable_type=">=", problem_type="-min")
sage: P.is_negative()
True

is_optimal(*x)

Check if given solution is feasible.

INPUT:

• anything that can be interpreted as a valid solution for this problem, i.e. a sequence of values for all decision variables

OUTPUT:

• True is the given solution is optimal, False otherwise

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (15, 5)
sage: P = InteractiveLPProblem(A, b, c, variable_type=">=")
sage: P.is_optimal(100, 200)
False
sage: P.is_optimal(500, 0)
True
sage: P.is_optimal(499, 3)
True
sage: P.is_optimal(501, -3)
False

is_primal()

Check if we consider this problem to be primal or dual.

This distinction affects only some automatically chosen variable names.

OUTPUT:

• boolean

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblem(A, b, c, ["C", "B"], variable_type=">=")
sage: P.is_primal()
True
sage: P.dual().is_primal()
False

m()

Return the number of constraints of self, i.e. $$m$$.

OUTPUT:

• an integer

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblem(A, b, c, ["C", "B"], variable_type=">=")
sage: P.n_constraints()
2
sage: P.m()
2

n()

Return the number of decision variables of self, i.e. $$n$$.

OUTPUT:

• an integer

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblem(A, b, c, ["C", "B"], variable_type=">=")
sage: P.n_variables()
2
sage: P.n()
2

n_constraints()

Return the number of constraints of self, i.e. $$m$$.

OUTPUT:

• an integer

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblem(A, b, c, ["C", "B"], variable_type=">=")
sage: P.n_constraints()
2
sage: P.m()
2

n_variables()

Return the number of decision variables of self, i.e. $$n$$.

OUTPUT:

• an integer

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblem(A, b, c, ["C", "B"], variable_type=">=")
sage: P.n_variables()
2
sage: P.n()
2

objective_coefficients()

Return coefficients of the objective of self, i.e. $$c$$.

OUTPUT:

• a vector

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblem(A, b, c, ["C", "B"], variable_type=">=")
sage: P.objective_coefficients()
(10, 5)
sage: P.c()
(10, 5)

objective_constant_term()

Return the constant term of the objective.

OUTPUT:

• a number

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblem(A, b, c, ["C", "B"], variable_type=">=")
sage: P.objective_constant_term()
0
sage: P.optimal_value()
6250
sage: P = InteractiveLPProblem(A, b, c, ["C", "B"],
....:       variable_type=">=", objective_constant_term=-1250)
sage: P.objective_constant_term()
-1250
sage: P.optimal_value()
5000

objective_value(*x)

Return the value of the objective on the given solution.

INPUT:

• anything that can be interpreted as a valid solution for this problem, i.e. a sequence of values for all decision variables

OUTPUT:

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblem(A, b, c, variable_type=">=")
sage: P.objective_value(100, 200)
2000

optimal_solution()

Return an optimal solution of self.

OUTPUT:

• a vector or None if there are no optimal solutions

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblem(A, b, c, ["C", "B"], variable_type=">=")
sage: P.optimal_solution()
(250, 750)

optimal_value()

Return the optimal value for self.

OUTPUT:

• a number if the problem is bounded, $$\pm \infty$$ if it is unbounded, or None if it is infeasible

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblem(A, b, c, ["C", "B"], variable_type=">=")
sage: P.optimal_value()
6250

plot(*args, **kwds)

Return a plot for solving self graphically.

INPUT:

• xmin, xmax, ymin, ymax – bounds for the axes, if not given, an attempt will be made to pick reasonable values
• alpha – (default: 0.2) determines how opaque are shadows

OUTPUT:

• a plot

This only works for problems with two decision variables. On the plot the black arrow indicates the direction of growth of the objective. The lines perpendicular to it are level curves of the objective. If there are optimal solutions, the arrow originates in one of them and the corresponding level curve is solid: all points of the feasible set on it are optimal solutions. Otherwise the arrow is placed in the center. If the problem is infeasible or the objective is zero, a plot of the feasible set only is returned.

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblem(A, b, c, ["C", "B"], variable_type=">=")
sage: p = P.plot()
sage: p.show()

In this case the plot works better with the following axes ranges:

sage: p = P.plot(0, 1000, 0, 1500)
sage: p.show()

plot_feasible_set(xmin=None, xmax=None, ymin=None, ymax=None, alpha=0.2)

Return a plot of the feasible set of self.

INPUT:

• xmin, xmax, ymin, ymax – bounds for the axes, if not given, an attempt will be made to pick reasonable values
• alpha – (default: 0.2) determines how opaque are shadows

OUTPUT:

• a plot

This only works for a problem with two decision variables. The plot shows boundaries of constraints with a shadow on one side for inequalities. If the feasible_set() is not empty and at least part of it is in the given boundaries, it will be shaded gray and $$F$$ will be placed in its middle.

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblem(A, b, c, ["C", "B"], variable_type=">=")
sage: p = P.plot_feasible_set()
sage: p.show()

In this case the plot works better with the following axes ranges:

sage: p = P.plot_feasible_set(0, 1000, 0, 1500)
sage: p.show()

problem_type()

Return the problem type.

Needs to be used together with is_negative.

OUTPUT:

• a string, one of "max", "min".

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblem(A, b, c, ["C", "B"], variable_type=">=")
sage: P.problem_type()
'max'
sage: P = InteractiveLPProblem(A, b, c, ["C", "B"], variable_type=">=", problem_type="-min")
sage: P.problem_type()
'min'

standard_form(transformation=False, **kwds)

Construct the LP problem in standard form equivalent to self.

INPUT:

• transformation – (default: False) if True, a map converting solutions of the problem in standard form to the original one will be returned as well
• you can pass (as keywords only) slack_variables, auxiliary_variable,objective_name to the constructor of InteractiveLPProblemStandardForm

OUTPUT:

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblem(A, b, c, variable_type=">=")
sage: DP = P.dual()
sage: DPSF = DP.standard_form()
sage: DPSF.b()
(-10, -5)
sage: DPSF.slack_variables()
(y3, y4)
sage: DPSF = DP.standard_form(slack_variables=["L", "F"])
sage: DPSF.slack_variables()
(L, F)
sage: DPSF, f = DP.standard_form(True)
sage: f
Vector space morphism represented by the matrix:
[1 0]
[0 1]
Domain: Vector space of dimension 2 over Rational Field
Codomain: Vector space of dimension 2 over Rational Field

A more complicated transformation map:

sage: P = InteractiveLPProblem(A, b, c, variable_type=["<=", ""],
....:                          objective_constant_term=42)
sage: PSF, f = P.standard_form(True)
sage: f
Vector space morphism represented by the matrix:
[-1  0]
[ 0  1]
[ 0 -1]
Domain: Vector space of dimension 3 over Rational Field
Codomain: Vector space of dimension 2 over Rational Field
sage: PSF.optimal_solution()
(0, 1000, 0)
sage: P.optimal_solution()
(0, 1000)
sage: P.is_optimal(PSF.optimal_solution())
Traceback (most recent call last):
...
TypeError: given input is not a solution for this problem
sage: P.is_optimal(f(PSF.optimal_solution()))
True
sage: PSF.optimal_value()
5042
sage: P.optimal_value()
5042

variable_types()

Return a tuple listing the variable types of all decision variables.

OUTPUT:

• a tuple of strings

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblem(A, b, c, ["C", "B"], variable_type=[">=", ""])
sage: P.variable_types()
('>=', '')

x()

Return decision variables of self, i.e. $$x$$.

OUTPUT:

• a vector

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblem(A, b, c, ["C", "B"], variable_type=">=")
sage: P.decision_variables()
(C, B)
sage: P.x()
(C, B)

class sage.numerical.interactive_simplex_method.InteractiveLPProblemStandardForm(A, b, c, x='x', problem_type='max', slack_variables=None, auxiliary_variable=None, base_ring=None, is_primal=True, objective_name=None, objective_constant_term=0)

Construct an LP (Linear Programming) problem in standard form.

Note

This class is for educational purposes only: if you want to solve Linear Programs efficiently, use MixedIntegerLinearProgram instead.

The used standard form is:

$\begin{split}\begin{array}{l} \pm \max cx \\ Ax \leq b \\ x \geq 0 \end{array}\end{split}$

INPUT:

• A – a matrix of constraint coefficients
• b – a vector of constraint constant terms
• c – a vector of objective coefficients
• x – (default: "x") a vector of decision variables or a string the base name giving
• problem_type – (default: "max") a string specifying the problem type: either "max" or "-max"
• slack_variables – (default: depends on style()) a vector of slack variables or a string giving the base name
• auxiliary_variable – (default: same as x parameter with adjoined "0" if it was given as a string, otherwise "x0") the auxiliary name, expected to be the same as the first decision variable for auxiliary problems
• base_ring – (default: the fraction field of a common ring for all input coefficients) a field to which all input coefficients will be converted
• is_primal – (default: True) whether this problem is primal or dual: each problem is of course dual to its own dual, this flag is mostly for internal use and affects default variable names only
• objective_name – a string or a symbolic expression for the objective used in dictionaries, default depends on style()
• objective_constant_term – (default: 0) a constant term of the objective

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)

Unlike InteractiveLPProblem, this class does not allow you to adjust types of constraints (they are always "<=") and variables (they are always ">="), and the problem type may only be "max" or "-max". You may give custom names to slack and auxiliary variables, but in most cases defaults should work:

sage: P.decision_variables()
(x1, x2)
sage: P.slack_variables()
(x3, x4)

add_constraint(coefficients, constant_term, slack_variable=None)

Return a new LP problem by adding a constraint toself.

INPUT:

• coefficients – coefficients of the new constraint
• constant_term – a constant term of the new constraint
• slack_variable – (default: depends on style()) a string giving the name of the slack variable of the new constraint

OUTPUT:

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: P.Abcx()
(
[1 1]
[3 1], (1000, 1500), (10, 5), (x1, x2)
)
sage: P.slack_variables()
(x3, x4)
sage: P1 = P.add_constraint(([2, 4]), 2000)
sage: P1.Abcx()
(
[1 1]
[3 1]
[2 4], (1000, 1500, 2000), (10, 5), (x1, x2)
)
sage: P1.slack_variables()
(x3, x4, x5)
sage: P2 = P.add_constraint(([2, 4]), 2000, slack_variable='c')
sage: P2.slack_variables()
(x3, x4, c)
sage: P3 = P.add_constraint(([2, 4, 6]), 2000)
Traceback (most recent call last):
...
TypeError: number of columns must be the same, not 2 and 3

auxiliary_problem(objective_name=None)

Construct the auxiliary problem for self.

INPUT:

• objective_name – a string or a symbolic expression for the objective used in dictionaries, default depends on style()

OUTPUT:

The auxiliary problem with the auxiliary variable $$x_0$$ is

$\begin{split}\begin{array}{l} \max - x_0 \\ - x_0 + A_i x \leq b_i \text{ for all } i \\ x \geq 0 \end{array}\ .\end{split}$

Such problems are used when the initial_dictionary() is infeasible.

EXAMPLES:

sage: A = ([1, 1], [3, 1], [-1, -1])
sage: b = (1000, 1500, -400)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: AP = P.auxiliary_problem()

auxiliary_variable()

Return the auxiliary variable of self.

Note that the auxiliary variable may or may not be among decision_variables().

OUTPUT:

EXAMPLES:

sage: A = ([1, 1], [3, 1], [-1, -1])
sage: b = (1000, 1500, -400)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: P.auxiliary_variable()
x0
sage: P.decision_variables()
(x1, x2)
sage: AP = P.auxiliary_problem()
sage: AP.auxiliary_variable()
x0
sage: AP.decision_variables()
(x0, x1, x2)

coordinate_ring()

Return the coordinate ring of self.

OUTPUT:

EXAMPLES:

sage: A = ([1, 1], [3, 1], [-1, -1])
sage: b = (1000, 1500, -400)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: P.coordinate_ring()
Multivariate Polynomial Ring in x0, x1, x2, x3, x4, x5
over Rational Field
sage: P.base_ring()
Rational Field
sage: P.auxiliary_variable()
x0
sage: P.decision_variables()
(x1, x2)
sage: P.slack_variables()
(x3, x4, x5)

dictionary(*x_B)

Construct a dictionary for self with given basic variables.

INPUT:

• basic variables for the dictionary to be constructed

OUTPUT:

Note

This is a synonym for self.revised_dictionary(x_B).dictionary(), but basic variables are mandatory.

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.dictionary("x1", "x2")
sage: D.basic_variables()
(x1, x2)

feasible_dictionary(auxiliary_dictionary)

Construct a feasible dictionary for self.

INPUT:

• auxiliary_dictionary – an optimal dictionary for the auxiliary_problem() of self with the optimal value $$0$$ and a non-basic auxiliary variable

OUTPUT:

EXAMPLES:

sage: A = ([1, 1], [3, 1], [-1, -1])
sage: b = (1000, 1500, -400)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: AP = P.auxiliary_problem()
sage: D = AP.initial_dictionary()
sage: D.enter(0)
sage: D.leave(5)
sage: D.update()
sage: D.enter(1)
sage: D.leave(0)
sage: D.update()
sage: D.is_optimal()
True
sage: D.objective_value()
0
sage: D.basic_solution()
(0, 400, 0)
sage: D = P.feasible_dictionary(D)
sage: D.is_optimal()
False
sage: D.is_feasible()
True
sage: D.objective_value()
4000
sage: D.basic_solution()
(400, 0)

final_dictionary()

Return the final dictionary of the simplex method applied to self.

See run_simplex_method() for the description of possibilities.

OUTPUT:

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.final_dictionary()
sage: D.is_optimal()
True

final_revised_dictionary()

Return the final dictionary of the revised simplex method applied to self.

See run_revised_simplex_method() for the description of possibilities.

OUTPUT:

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.final_revised_dictionary()
sage: D.is_optimal()
True

initial_dictionary()

Construct the initial dictionary of self.

The initial dictionary “defines” slack_variables() in terms of the decision_variables(), i.e. it has slack variables as basic ones.

OUTPUT:

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.initial_dictionary()

inject_variables(scope=None, verbose=True)

Inject variables of self into scope.

INPUT:

• scope – namespace (default: global)
• verbose – if True (default), names of injected variables will be printed

OUTPUT:

• none

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: P.inject_variables()
Defining x0, x1, x2, x3, x4
sage: 3*x1 + x2
x2 + 3*x1

objective_name()

Return the objective name used in dictionaries for this problem.

OUTPUT:

• a symbolic expression

EXAMPLES:

sage: A = ([1, 1], [3, 1], [-1, -1])
sage: b = (1000, 1500, -400)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: P.objective_name()
z
sage: sage.numerical.interactive_simplex_method.style("Vanderbei")
'Vanderbei'
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: P.objective_name()
zeta
sage: sage.numerical.interactive_simplex_method.style("UAlberta")
'UAlberta'
sage: P = InteractiveLPProblemStandardForm(A, b, c, objective_name="custom")
sage: P.objective_name()
custom

revised_dictionary(*x_B)

Construct a revised dictionary for self.

INPUT:

OUTPUT:

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.revised_dictionary("x1", "x2")
sage: D.basic_variables()
(x1, x2)

If basic variables are not given the initial dictionary is constructed:

sage: P.revised_dictionary().basic_variables()
(x3, x4)
sage: P.initial_dictionary().basic_variables()
(x3, x4)

Unless it is infeasible, in which case a feasible dictionary for the auxiliary problem is constructed:

sage: A = ([1, 1], [3, 1], [-1,-1])
sage: b = (1000, 1500, -400)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: P.initial_dictionary().is_feasible()
False
sage: P.revised_dictionary().basic_variables()
(x3, x4, x0)

run_revised_simplex_method()

Apply the revised simplex method and return all steps.

OUTPUT:

Note

You can access the final_revised_dictionary(), which can be one of the following:

• an optimal dictionary with the auxiliary_variable() among basic_variables() and a non-zero optimal value indicating that self is infeasible;
• a non-optimal dictionary that has marked entering variable for which there is no choice of the leaving variable, indicating that self is unbounded;
• an optimal dictionary.

EXAMPLES:

sage: A = ([1, 1], [3, 1], [-1, -1])
sage: b = (1000, 1500, -400)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: P.run_revised_simplex_method()
\begin{equation*}
...
\end{equation*}
Entering: $x_{1}$. Leaving: $x_{0}$.
\begin{equation*}
...
\end{equation*}
Entering: $x_{5}$. Leaving: $x_{4}$.
\begin{equation*}
...
\end{equation*}
Entering: $x_{2}$. Leaving: $x_{3}$.
\begin{equation*}
...
\end{equation*}
The optimal value: $6250$. An optimal solution: $\left(250,\,750\right)$.

run_simplex_method()

Apply the simplex method and return all steps and intermediate states.

OUTPUT:

Note

You can access the final_dictionary(), which can be one of the following:

• an optimal dictionary for the auxiliary_problem() with a non-zero optimal value indicating that self is infeasible;
• a non-optimal dictionary for self that has marked entering variable for which there is no choice of the leaving variable, indicating that self is unbounded;
• an optimal dictionary for self.

EXAMPLES:

sage: A = ([1, 1], [3, 1], [-1, -1])
sage: b = (1000, 1500, -400)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: P.run_simplex_method()
\begin{equation*}
...
\end{equation*}
The initial dictionary is infeasible, solving auxiliary problem.
...
Entering: $x_{0}$. Leaving: $x_{5}$.
...
Entering: $x_{1}$. Leaving: $x_{0}$.
...
Back to the original problem.
...
Entering: $x_{5}$. Leaving: $x_{4}$.
...
Entering: $x_{2}$. Leaving: $x_{3}$.
...
The optimal value: $6250$. An optimal solution: $\left(250,\,750\right)$.

slack_variables()

Return slack variables of self.

Slack variables are differences between the constant terms and left hand sides of the constraints.

If you want to give custom names to slack variables, you have to do so during construction of the problem.

OUTPUT:

• a tuple

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: P.slack_variables()
(x3, x4)
sage: P = InteractiveLPProblemStandardForm(A, b, c, ["C", "B"],
....:     slack_variables=["L", "F"])
sage: P.slack_variables()
(L, F)

class sage.numerical.interactive_simplex_method.LPAbstractDictionary

Abstract base class for dictionaries for LP problems.

Instantiating this class directly is meaningless, see LPDictionary and LPRevisedDictionary for useful extensions.

add_row(nonbasic_coefficients, constant, basic_variable=None)

Return a dictionary with an additional row based on a given dictionary.

INPUT:

• nonbasic_coefficients– a list of the coefficients for the new row (with which nonbasic variables are subtracted in the relation for the new basic variable)
• constant– the constant term for the new row
• basic_variable– (default: depends on style()) a string giving the name of the basic variable of the new row

OUTPUT:

• a new dictionary of the same class

EXAMPLES:

sage: A = ([-1, 1, 7], [8, 2, 13], [34, 17, 12])
sage: b = (2, 17, 6)
sage: c = (55/10, 21/10, 14/30)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.dictionary("x1", "x2", "x4")
sage: D1 = D.add_row([7, 11, 19], 42, basic_variable='c')
sage: D1.row_coefficients("c")
(7, 11, 19)

base_ring()

Return the base ring of self, i.e. the ring of coefficients.

OUTPUT:

• a ring

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.initial_dictionary()
sage: D.base_ring()
Rational Field
sage: D = P.revised_dictionary()
sage: D.base_ring()
Rational Field

basic_solution(include_slack_variables=False)

Return the basic solution of self.

The basic solution associated to a dictionary is obtained by setting to zero all nonbasic_variables(), in which case basic_variables() have to be equal to constant_terms() in equations. It may refer to values of decision_variables() only or include slack_variables() as well.

INPUT:

• include_slack_variables – (default: False) if True, values of slack variables will be appended at the end

OUTPUT:

• a vector

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.initial_dictionary()
sage: D.basic_solution()
(0, 0)
sage: D.basic_solution(True)
(0, 0, 1000, 1500)
sage: D = P.revised_dictionary()
sage: D.basic_solution()
(0, 0)
sage: D.basic_solution(True)
(0, 0, 1000, 1500)

basic_variables()

Return the basic variables of self.

OUTPUT:

• a vector

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.initial_dictionary()
sage: D.basic_variables()
(x3, x4)

column_coefficients(v)

Return the coefficients of a nonbasic variable.

INPUT:

• v – a nonbasic variable of self, can be given as a string, an actual variable, or an integer interpreted as the index of a variable

OUTPUT:

• a vector of coefficients of a nonbasic variable

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.revised_dictionary()
sage: D.column_coefficients(1)
(1, 3)

constant_terms()

Return the constant terms of relations of self.

OUTPUT:

• a vector.

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.initial_dictionary()
sage: D.constant_terms()
(1000, 1500)

coordinate_ring()

Return the coordinate ring of self.

OUTPUT:

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.initial_dictionary()
sage: D.coordinate_ring()
Multivariate Polynomial Ring in x0, x1, x2, x3, x4
over Rational Field
sage: D = P.revised_dictionary()
sage: D.coordinate_ring()
Multivariate Polynomial Ring in x0, x1, x2, x3, x4
over Rational Field

dual_ratios()

Return ratios used to determine the entering variable based on leaving.

OUTPUT:

• A list of pairs $$(r_j, x_j)$$ where $$x_j$$ is a non-basic variable and $$r_j = c_j / a_{ij}$$ is the ratio of the objective coefficient $$c_j$$ to the coefficient $$a_{ij}$$ of $$x_j$$ in the relation for the leaving variable $$x_i$$:

$x_i = b_i - \cdots - a_{ij} x_j - \cdots.$

The order of pairs matches the order of nonbasic_variables(), but only $$x_j$$ with negative $$a_{ij}$$ are considered.

EXAMPLES:

sage: A = ([1, 1], [3, 1], [-1, -1])
sage: b = (1000, 1500, -400)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.dictionary(2, 3, 5)
sage: D.leave(3)
sage: D.dual_ratios()
[(5/2, x1), (5, x4)]
sage: D = P.revised_dictionary(2, 3, 5)
sage: D.leave(3)
sage: D.dual_ratios()
[(5/2, x1), (5, x4)]

enter(v)

Set v as the entering variable of self.

INPUT:

• v – a non-basic variable of self, can be given as a string, an actual variable, or an integer interpreted as the index of a variable. It is also possible to enter None to reset choice.

OUTPUT:

• none, but the selected variable will be used as entering by methods that require an entering variable and the corresponding column will be typeset in green

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.initial_dictionary()
sage: D.enter("x1")

We can also use indices of variables:

sage: D.enter(1)

Or variable names without quotes after injecting them:

sage: P.inject_variables()
Defining x0, x1, x2, x3, x4
sage: D.enter(x1)

The same works for revised dictionaries as well:

sage: D = P.revised_dictionary()
sage: D.enter(x1)

entering()

Return the currently chosen entering variable.

OUTPUT:

• a variable if the entering one was chosen, otherwise None

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.initial_dictionary()
sage: D.entering() is None
True
sage: D.enter(1)
sage: D.entering()
x1

entering_coefficients()

Return coefficients of the entering variable.

OUTPUT:

• a vector

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.initial_dictionary()
sage: D.enter(1)
sage: D.entering_coefficients()
(1, 3)

is_dual_feasible()

Check if self is dual feasible.

OUTPUT:

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.initial_dictionary()
sage: D.is_dual_feasible()
False
sage: D = P.revised_dictionary()
sage: D.is_dual_feasible()
False

is_feasible()

Check if self is feasible.

OUTPUT:

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.initial_dictionary()
sage: D.is_feasible()
True
sage: D = P.revised_dictionary()
sage: D.is_feasible()
True

is_optimal()

Check if self is optimal.

OUTPUT:

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.initial_dictionary()
sage: D.is_optimal()
False
sage: D = P.revised_dictionary()
sage: D.is_optimal()
False
sage: D = P.revised_dictionary(1, 2)
sage: D.is_optimal()
True

leave(v)

Set v as the leaving variable of self.

INPUT:

• v – a basic variable of self, can be given as a string, an actual variable, or an integer interpreted as the index of a variable. It is also possible to leave None to reset choice.

OUTPUT:

• none, but the selected variable will be used as leaving by methods that require a leaving variable and the corresponding row will be typeset in red

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.initial_dictionary()
sage: D.leave("x4")

We can also use indices of variables:

sage: D.leave(4)

Or variable names without quotes after injecting them:

sage: P.inject_variables()
Defining x0, x1, x2, x3, x4
sage: D.leave(x4)

The same works for revised dictionaries as well:

sage: D = P.revised_dictionary()
sage: D.leave(x4)

leaving()

Return the currently chosen leaving variable.

OUTPUT:

• a variable if the leaving one was chosen, otherwise None

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.initial_dictionary()
sage: D.leaving() is None
True
sage: D.leave(4)
sage: D.leaving()
x4

leaving_coefficients()

Return coefficients of the leaving variable.

OUTPUT:

• a vector

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.dictionary(2, 3)
sage: D.leave(3)
sage: D.leaving_coefficients()
(-2, -1)

The same works for revised dictionaries as well:

sage: D = P.revised_dictionary(2, 3)
sage: D.leave(3)
sage: D.leaving_coefficients()
(-2, -1)

nonbasic_variables()

Return non-basic variables of self.

OUTPUT:

• a vector

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.initial_dictionary()
sage: D.nonbasic_variables()
(x1, x2)

objective_coefficients()

Return coefficients of the objective of self.

OUTPUT:

• a vector

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.initial_dictionary()
sage: D.objective_coefficients()
(10, 5)

objective_name()

Return the objective name of self.

OUTPUT:

• a symbolic expression

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.initial_dictionary()
sage: D.objective_name()
z

objective_value()

Return the value of the objective at the basic_solution() of self.

OUTPUT:

• a number

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.initial_dictionary()
sage: D.objective_value()
0

possible_dual_simplex_method_steps()

Return possible dual simplex method steps for self.

OUTPUT:

• A list of pairs (leaving, entering), where leaving is a basic variable that may leave() and entering is a list of non-basic variables that may enter() when leaving leaves. Note that entering may be empty, indicating that the problem is infeasible (since the dual one is unbounded).

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.dictionary(2, 3)
sage: D.possible_dual_simplex_method_steps()
[(x3, [x1])]
sage: D = P.revised_dictionary(2, 3)
sage: D.possible_dual_simplex_method_steps()
[(x3, [x1])]

possible_entering()

Return possible entering variables for self.

OUTPUT:

• a list of non-basic variables of self that can enter() on the next step of the (dual) simplex method

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.initial_dictionary()
sage: D.possible_entering()
[x1, x2]
sage: D = P.revised_dictionary()
sage: D.possible_entering()
[x1, x2]

possible_leaving()

Return possible leaving variables for self.

OUTPUT:

• a list of basic variables of self that can leave() on the next step of the (dual) simplex method

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.initial_dictionary()
sage: D.enter(1)
sage: D.possible_leaving()
[x4]
sage: D = P.revised_dictionary()
sage: D.enter(1)
sage: D.possible_leaving()
[x4]

possible_simplex_method_steps()

Return possible simplex method steps for self.

OUTPUT:

• A list of pairs (entering, leaving), where entering is a non-basic variable that may enter() and leaving is a list of basic variables that may leave() when entering enters. Note that leaving may be empty, indicating that the problem is unbounded.

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.initial_dictionary()
sage: D.possible_simplex_method_steps()
[(x1, [x4]), (x2, [x3])]
sage: D = P.revised_dictionary()
sage: D.possible_simplex_method_steps()
[(x1, [x4]), (x2, [x3])]

ratios()

Return ratios used to determine the leaving variable based on entering.

OUTPUT:

• A list of pairs $$(r_i, x_i)$$ where $$x_i$$ is a basic variable and $$r_i = b_i / a_{ik}$$ is the ratio of the constant term $$b_i$$ to the coefficient $$a_{ik}$$ of the entering variable $$x_k$$ in the relation for $$x_i$$:

$x_i = b_i - \cdots - a_{ik} x_k - \cdots.$

The order of pairs matches the order of basic_variables(), but only $$x_i$$ with positive $$a_{ik}$$ are considered.

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.initial_dictionary()
sage: D.enter(1)
sage: D.ratios()
[(1000, x3), (500, x4)]
sage: D = P.revised_dictionary()
sage: D.enter(1)
sage: D.ratios()
[(1000, x3), (500, x4)]

row_coefficients(v)

Return the coefficients of the basic variable v.

These are the coefficients with which nonbasic variables are subtracted in the relation for v.

INPUT:

• v – a basic variable of self, can be given as a string, an actual variable, or an integer interpreted as the index of a variable

OUTPUT:

• a vector of coefficients of a basic variable

EXAMPLES:

sage: A = ([-1, 1], [8, 2])
sage: b = (2, 17)
sage: c = (55/10, 21/10)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.final_dictionary()
sage: D.row_coefficients("x1")
(1/10, -1/5)

We can also use indices of variables:

sage: D.row_coefficients(1)
(1/10, -1/5)

Or use variable names without quotes after injecting them:

sage: P.inject_variables()
Defining x0, x1, x2, x3, x4
sage: D.row_coefficients(x1)
(1/10, -1/5)

run_dual_simplex_method()

Apply the dual simplex method and return all steps/intermediate states.

If either entering or leaving variables were already set, they will be used.

OUTPUT:

EXAMPLES:

sage: A = ([1, 1], [3, 1], [-1, -1])
sage: b = (1000, 1500, -400)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.initial_dictionary()
sage: D.run_dual_simplex_method()
Traceback (most recent call last):
...
ValueError: leaving variables can be determined for feasible
dictionaries with a set entering variable or for dual feasible
dictionaries

sage: D = P.dictionary(2, 3, 5)
sage: D.is_dual_feasible()
True
sage: D.is_optimal()
False
sage: D.run_dual_simplex_method()
\begin{equation*}
...
\end{equation*}
Leaving: $x_{3}$. Entering: $x_{1}$.
\begin{equation*}
...
\end{equation*}
sage: D.is_optimal()
True

This method detects infeasible problems:

sage: A = ([1, 0],)
sage: b = (-1,)
sage: c = (0, -1)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.initial_dictionary()
sage: D.run_dual_simplex_method()
\begin{equation*}
...
\end{equation*}
The problem is infeasible because of $x_{3}$ constraint.

run_simplex_method()

Apply the simplex method and return all steps and intermediate states.

If either entering or leaving variables were already set, they will be used.

OUTPUT:

EXAMPLES:

sage: A = ([1, 1], [3, 1], [-1, -1])
sage: b = (1000, 1500, -400)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.initial_dictionary()
sage: D.run_simplex_method()
Traceback (most recent call last):
...
ValueError: entering variables can be determined for feasible
dictionaries or for dual feasible dictionaries with a set leaving
variable

sage: D = P.dictionary(1, 3, 4)
sage: D.is_feasible()
True
sage: D.is_optimal()
False
sage: D.run_simplex_method()
\begin{equation*}
...
\end{equation*}
Entering: $x_{5}$. Leaving: $x_{4}$.
\begin{equation*}
...
\end{equation*}
Entering: $x_{2}$. Leaving: $x_{3}$.
\begin{equation*}
...
\end{equation*}
sage: D.is_optimal()
True

This method detects unbounded problems:

sage: A = ([1, 0],)
sage: b = (1,)
sage: c = (0, 1)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.initial_dictionary()
sage: D.run_simplex_method()
\begin{equation*}
...
\end{equation*}
The problem is unbounded in $x_{2}$ direction.

update()

Update self using previously set entering and leaving variables.

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.initial_dictionary()
sage: D.objective_value()
0
sage: D.enter("x1")
sage: D.leave("x4")
sage: D.update()
sage: D.objective_value()
5000

class sage.numerical.interactive_simplex_method.LPDictionary(A, b, c, objective_value, basic_variables, nonbasic_variables, objective_name)

Construct a dictionary for an LP problem.

A dictionary consists of the following data:

$\begin{split}\begin{array}{|l|} \hline x_B = b - A x_N\\ \hline z = z^* + c x_N\\ \hline \end{array}\end{split}$

INPUT:

• A – a matrix of relation coefficients
• b – a vector of relation constant terms
• c – a vector of objective coefficients
• objective_value – current value of the objective $$z^*$$
• basic_variables – a list of basic variables $$x_B$$
• nonbasic_variables – a list of non-basic variables $$x_N$$
• objective_name – a “name” for the objective $$z$$

OUTPUT:

Note

This constructor does not check correctness of input, as it is intended to be used internally by InteractiveLPProblemStandardForm.

EXAMPLES:

The intended way to use this class is indirect:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.initial_dictionary()
sage: D
LP problem dictionary (use typeset mode to see details)

But if you want you can create a dictionary without starting with an LP problem, here is construction of the same dictionary as above:

sage: A = matrix(QQ, ([1, 1], [3, 1]))
sage: b = vector(QQ, (1000, 1500))
sage: c = vector(QQ, (10, 5))
sage: R = PolynomialRing(QQ, "x1, x2, x3, x4", order="neglex")
sage: from sage.numerical.interactive_simplex_method \
....:     import LPDictionary
sage: D2 = LPDictionary(A, b, c, 0, R.gens()[2:], R.gens()[:2], "z")
sage: D2 == D
True

ELLUL(entering, leaving)

Perform the Enter-Leave-LaTeX-Update-LaTeX step sequence on self.

INPUT:

• entering – the entering variable
• leaving – the leaving variable

OUTPUT:

• a string with LaTeX code for self before and after update

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.initial_dictionary()
sage: D.ELLUL("x1", "x4")
doctest:...: DeprecationWarning: ELLUL is deprecated, please use separate enter-leave-update and output commands
See http://trac.sagemath.org/19097 for details.
\renewcommand{\arraystretch}{1.5} %notruncate
\begin{array}{|rcrcrcr|}
\hline
x_{3} \mspace{-6mu}&\mspace{-6mu} = \mspace{-6mu}&\mspace{-6mu} 1000 \mspace{-6mu}&\mspace{-6mu} - \mspace{-6mu}&\color{green}\mspace{-6mu} x_{1} \mspace{-6mu}&\mspace{-6mu} - \mspace{-6mu}&\mspace{-6mu} x_{2}\\
\color{red}x_{4} \mspace{-6mu}&\color{red}\mspace{-6mu} = \mspace{-6mu}&\color{red}\mspace{-6mu} 1500 \mspace{-6mu}&\color{red}\mspace{-6mu} - \mspace{-6mu}&\color{blue}\mspace{-6mu} 3 x_{1} \mspace{-6mu}&\color{red}\mspace{-6mu} - \mspace{-6mu}&\color{red}\mspace{-6mu} x_{2}\\
\hline
z \mspace{-6mu}&\mspace{-6mu} = \mspace{-6mu}&\mspace{-6mu} 0 \mspace{-6mu}&\mspace{-6mu} + \mspace{-6mu}&\color{green}\mspace{-6mu} 10 x_{1} \mspace{-6mu}&\mspace{-6mu} + \mspace{-6mu}&\mspace{-6mu} 5 x_{2}\\
\hline
\\
\hline
x_{3} \mspace{-6mu}&\mspace{-6mu} = \mspace{-6mu}&\mspace{-6mu} 500 \mspace{-6mu}&\mspace{-6mu} + \mspace{-6mu}&\mspace{-6mu} \frac{1}{3} x_{4} \mspace{-6mu}&\mspace{-6mu} - \mspace{-6mu}&\mspace{-6mu} \frac{2}{3} x_{2}\\
x_{1} \mspace{-6mu}&\mspace{-6mu} = \mspace{-6mu}&\mspace{-6mu} 500 \mspace{-6mu}&\mspace{-6mu} - \mspace{-6mu}&\mspace{-6mu} \frac{1}{3} x_{4} \mspace{-6mu}&\mspace{-6mu} - \mspace{-6mu}&\mspace{-6mu} \frac{1}{3} x_{2}\\
\hline
z \mspace{-6mu}&\mspace{-6mu} = \mspace{-6mu}&\mspace{-6mu} 5000 \mspace{-6mu}&\mspace{-6mu} - \mspace{-6mu}&\mspace{-6mu} \frac{10}{3} x_{4} \mspace{-6mu}&\mspace{-6mu} + \mspace{-6mu}&\mspace{-6mu} \frac{5}{3} x_{2}\\
\hline
\end{array}

This is how the above output looks when rendered:

$\begin{split}\renewcommand{\arraystretch}{1.5} \begin{array}{|rcrcrcr|} \hline x_{3} \mspace{-6mu}&\mspace{-6mu} = \mspace{-6mu}&\mspace{-6mu} 1000 \mspace{-6mu}&\mspace{-6mu} - \mspace{-6mu}&\color{green}\mspace{-6mu} x_{1} \mspace{-6mu}&\mspace{-6mu} - \mspace{-6mu}&\mspace{-6mu} x_{2}\\ \color{red}x_{4} \mspace{-6mu}&\color{red}\mspace{-6mu} = \mspace{-6mu}&\color{red}\mspace{-6mu} 1500 \mspace{-6mu}&\color{red}\mspace{-6mu} - \mspace{-6mu}&\color{blue}\mspace{-6mu} 3 x_{1} \mspace{-6mu}&\color{red}\mspace{-6mu} - \mspace{-6mu}&\color{red}\mspace{-6mu} x_{2}\\ \hline z \mspace{-6mu}&\mspace{-6mu} = \mspace{-6mu}&\mspace{-6mu} 0 \mspace{-6mu}&\mspace{-6mu} + \mspace{-6mu}&\color{green}\mspace{-6mu} 10 x_{1} \mspace{-6mu}&\mspace{-6mu} + \mspace{-6mu}&\mspace{-6mu} 5 x_{2}\\ \hline \\ \hline x_{3} \mspace{-6mu}&\mspace{-6mu} = \mspace{-6mu}&\mspace{-6mu} 500 \mspace{-6mu}&\mspace{-6mu} + \mspace{-6mu}&\mspace{-6mu} \frac{1}{3} x_{4} \mspace{-6mu}&\mspace{-6mu} - \mspace{-6mu}&\mspace{-6mu} \frac{2}{3} x_{2}\\ x_{1} \mspace{-6mu}&\mspace{-6mu} = \mspace{-6mu}&\mspace{-6mu} 500 \mspace{-6mu}&\mspace{-6mu} - \mspace{-6mu}&\mspace{-6mu} \frac{1}{3} x_{4} \mspace{-6mu}&\mspace{-6mu} - \mspace{-6mu}&\mspace{-6mu} \frac{1}{3} x_{2}\\ \hline z \mspace{-6mu}&\mspace{-6mu} = \mspace{-6mu}&\mspace{-6mu} 5000 \mspace{-6mu}&\mspace{-6mu} - \mspace{-6mu}&\mspace{-6mu} \frac{10}{3} x_{4} \mspace{-6mu}&\mspace{-6mu} + \mspace{-6mu}&\mspace{-6mu} \frac{5}{3} x_{2}\\ \hline \end{array}\end{split}$

The column of the entering variable is green, while the row of the leaving variable is red in the original dictionary state on the top. The new state after the update step is shown on the bottom.

add_row(nonbasic_coefficients, constant, basic_variable=None)

Return a dictionary with an additional row based on a given dictionary.

INPUT:

• nonbasic_coefficients– a list of the coefficients for the new row (with which nonbasic variables are subtracted in the relation for the new basic variable)
• constant– the constant term for the new row
• basic_variable– (default: depends on style()) a string giving the name of the basic variable of the new row

OUTPUT:

EXAMPLES:

sage: A = ([-1, 1, 7], [8, 2, 13], [34, 17, 12])
sage: b = (2, 17, 6)
sage: c = (55/10, 21/10, 14/30)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.dictionary("x1", "x2", "x4")
sage: D1 = D.add_row([7, 11, 19], 42, basic_variable='c')
sage: D1.row_coefficients("c")
(7, 11, 19)
sage: D1.constant_terms()[-1]
42
sage: D1.basic_variables()[-1]
c

basic_variables()

Return the basic variables of self.

OUTPUT:

• a vector

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.initial_dictionary()
sage: D.basic_variables()
(x3, x4)

column_coefficients(v)

Return coefficients of a nonbasic variable.

INPUT:

• v – a nonbasic variable of self, can be given as a string, an actual variable, or an integer interpreted as the index of a variable

OUTPUT:

• a vector

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.initial_dictionary()
sage: D.column_coefficients(1)
(1, 3)

constant_terms()

Return the constant terms of relations of self.

OUTPUT:

• a vector.

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.initial_dictionary()
sage: D.constant_terms()
(1000, 1500)

nonbasic_variables()

Return non-basic variables of self.

OUTPUT:

• a vector

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.initial_dictionary()
sage: D.nonbasic_variables()
(x1, x2)

objective_coefficients()

Return coefficients of the objective of self.

OUTPUT:

• a vector

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.initial_dictionary()
sage: D.objective_coefficients()
(10, 5)

objective_name()

Return the objective name of self.

OUTPUT:

• a symbolic expression

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.initial_dictionary()
sage: D.objective_name()
z

objective_value()

Return the value of the objective at the basic_solution() of self.

OUTPUT:

• a number

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.initial_dictionary()
sage: D.objective_value()
0

row_coefficients(v)

Return the coefficients of the basic variable v.

These are the coefficients with which nonbasic variables are subtracted in the relation for v.

INPUT:

• v – a basic variable of self, can be given as a string, an actual variable, or an integer interpreted as the index of a variable

OUTPUT:

• a vector of coefficients of a basic variable

EXAMPLES:

sage: A = ([-1, 1], [8, 2])
sage: b = (2, 17)
sage: c = (55/10, 21/10)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.final_dictionary()
sage: D.row_coefficients("x1")
(1/10, -1/5)

We can also use indices of variables:

sage: D.row_coefficients(1)
(1/10, -1/5)

Or use variable names without quotes after injecting them:

sage: P.inject_variables()
Defining x0, x1, x2, x3, x4
sage: D.row_coefficients(x1)
(1/10, -1/5)

update()

Update self using previously set entering and leaving variables.

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.initial_dictionary()
sage: D.objective_value()
0
sage: D.enter("x1")
sage: D.leave("x4")
sage: D.update()
sage: D.objective_value()
5000

class sage.numerical.interactive_simplex_method.LPRevisedDictionary(problem, basic_variables)

Construct a revised dictionary for an LP problem.

INPUT:

OUTPUT:

A revised dictionary encodes the same relations as a regular dictionary, but stores only what is “necessary to efficiently compute data for the simplex method”.

Let the original problem be

$\begin{split}\begin{array}{l} \pm \max cx \\ Ax \leq b \\ x \geq 0 \end{array}\end{split}$

Let $$\bar{x}$$ be the vector of decision_variables() $$x$$ followed by the slack_variables(). Let $$\bar{c}$$ be the vector of objective_coefficients() $$c$$ followed by zeroes for all slack variables. Let $$\bar{A} = (A | I)$$ be the matrix of constraint_coefficients() $$A$$ augmented by the identity matrix as columns corresponding to the slack variables. Then the problem above can be written as

$\begin{split}\begin{array}{l} \pm \max \bar{c} \bar{x} \\ \bar{A} \bar{x} = b \\ \bar{x} \geq 0 \end{array}\end{split}$

and any dictionary is a system of equations equivalent to $$\bar{A} \bar{x} = b$$, but resolved for basic_variables() $$x_B$$ in terms of nonbasic_variables() $$x_N$$ together with the expression for the objective in terms of $$x_N$$. Let c_B() and c_N() be vectors “splitting $$\bar{c}$$ into basic and non-basic parts”. Let B() and A_N() be the splitting of $$\bar{A}$$. Then the corresponding dictionary is

$\begin{split}\begin{array}{|l|} \hline x_B = B^{-1} b - B^{-1} A_N x_N\\ \hline z = y b + \left(c_N - y^T A_N\right) x_N\\ \hline \end{array}\end{split}$

where $$y = c_B^T B^{-1}$$. To proceed with the simplex method, it is not necessary to compute all entries of this dictionary. On the other hand, any entry is easy to compute, if you know $$B^{-1}$$, so we keep track of it through the update steps.

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: from sage.numerical.interactive_simplex_method \
....:     import LPRevisedDictionary
sage: D = LPRevisedDictionary(P, [1, 2])
sage: D.basic_variables()
(x1, x2)
sage: D
LP problem dictionary (use typeset mode to see details)

The same dictionary can be constructed through the problem:

sage: P.revised_dictionary(1, 2) == D
True

When this dictionary is typeset, you will see two tables like these ones:

$\begin{split}\renewcommand{\arraystretch}{1.500000} \begin{array}{l} \begin{array}{l|r|rr||r||r} x_B & c_B & & \mspace{-16mu} B^{-1} & y & B^{-1} b \\ \hline x_{1} & 10 & -\frac{1}{2} & \frac{1}{2} & \frac{5}{2} & 250 \\ x_{2} & 5 & \frac{3}{2} & -\frac{1}{2} & \frac{5}{2} & 750 \\ \end{array}\\ \\ \begin{array}{r|rr} x_N & x_{3} & x_{4} \\ \hline c_N^T & 0 & 0 \\ \hline y^T A_N & \frac{5}{2} & \frac{5}{2} \\ \hline c_N^T - y^T A_N & -\frac{5}{2} & -\frac{5}{2} \\ \end{array} \end{array}\end{split}$

More details will be shown if entering and leaving variables are set, but in any case the top table shows $$B^{-1}$$ and a few extra columns, while the bottom one shows several rows: these are related to columns and rows of dictionary entries.

A(v)

Return the column of constraint coefficients corresponding to v.

INPUT:

• v – a variable, its name, or its index

OUTPUT:

• a vector

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.revised_dictionary()
sage: D.A(1)
(1, 3)
sage: D.A(0)
(-1, -1)
sage: D.A("x3")
(1, 0)

A_N()

Return the $$A_N$$ matrix, constraint coefficients of non-basic variables.

OUTPUT:

• a matrix

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.revised_dictionary()
sage: D.A_N()
[1 1]
[3 1]

B()

Return the $$B$$ matrix, i.e. constraint coefficients of basic variables.

OUTPUT:

• a matrix

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.revised_dictionary(1, 2)
sage: D.B()
[1 1]
[3 1]

B_inverse()

Return the inverse of the B() matrix.

This inverse matrix is stored and computed during dictionary update in a more efficient way than generic inversion.

OUTPUT:

• a matrix

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.revised_dictionary(1, 2)
sage: D.B_inverse()
[-1/2  1/2]
[ 3/2 -1/2]

E()

Return the eta matrix between self and the next dictionary.

OUTPUT:

• a matrix

If $$B_{\mathrm{old}}$$ is the current matrix $$B$$ and $$B_{\mathrm{new}}$$ is the $$B$$ matrix of the next dictionary (after the update step), then $$B_{\mathrm{new}} = B_{\mathrm{old}} E$$.

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.revised_dictionary()
sage: D.enter(1)
sage: D.leave(4)
sage: D.E()
[1 1]
[0 3]

E_inverse()

Return the inverse of the matrix E().

This inverse matrix is computed in a more efficient way than generic inversion.

OUTPUT:

• a matrix

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.revised_dictionary()
sage: D.enter(1)
sage: D.leave(4)
sage: D.E_inverse()
[   1 -1/3]
[   0  1/3]

add_row(nonbasic_coefficients, constant, basic_variable=None)

Return a dictionary with an additional row based on a given dictionary.

The implementation of this method for revised dictionaries adds a new inequality constraint to the problem, in which the given $$basic_variable$$ becomes the slack variable. The resulting dictionary (with $$basic_variable$$ added to the basis) will have the given $$nonbasic_coefficients$$ and $$constant$$ as a new row.

INPUT:

• nonbasic_coefficients– a list of the coefficients for the new row (with which nonbasic variables are subtracted in the relation for the new basic variable)
• constant– the constant term for the new row
• basic_variable– (default: depends on style()) a string giving the name of the basic variable of the new row

OUTPUT:

EXAMPLES:

sage: A = ([-1, 1111, 3, 17], [8, 222, 7, 6],
....: [3, 7, 17, 5], [9, 5, 7, 3])
sage: b = (2, 17, 11, 27)
sage: c = (5/133, 1/10, 1/18, 47/3)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.final_revised_dictionary()
sage: D1 = D.add_row([7, 11, 13, 9], 42)
sage: D1.row_coefficients("x9")
(7, 11, 13, 9)
sage: D1.constant_terms()[-1]
42
sage: D1.basic_variables()[-1]
x9

sage: A = ([-9, 7, 48, 31, 23], [5, 2, 9, 13, 98],
....: [14, 15, 97, 49, 1], [9, 5, 7, 3, 17],
....: [119, 7, 121, 5, 111])
sage: b = (33, 27, 1, 272, 61)
sage: c = (51/133, 1/100, 149/18, 47/37, 13/17)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.revised_dictionary("x1", "x2", "x3", "x4", "x5")
sage: D2 = D.add_row([5 ,7, 11, 13, 9], 99, basic_variable='c')
sage: D2.row_coefficients("c")
(5, 7, 11, 13, 9)
sage: D2.constant_terms()[-1]
99
sage: D2.basic_variables()[-1]
c

sage: D = P.revised_dictionary(0, 1, 2, 3, 4)
sage: D.add_row([1, 2, 3, 4, 5, 6], 0)
Traceback (most recent call last):
...
ValueError: the sum of coefficients of nonbasic slack variables has
to be equal to -1 when inserting a row into a dictionary for the
auxiliary problem
sage: D3 = D.add_row([1, 2, 3, 4, 5, -15], 0)
sage: D3.row_coefficients(11)
(1, 2, 3, 4, 5, -15)

basic_indices()

Return the basic indices of self.

Note

Basic indices are indices of basic_variables() in the list of generators of the coordinate_ring() of the problem() of self, they may not coincide with the indices of variables which are parts of their names. (They will for the default indexed names.)

OUTPUT:

• a list.

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.revised_dictionary()
sage: D.basic_indices()
[3, 4]

basic_variables()

Return the basic variables of self.

OUTPUT:

• a vector

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.revised_dictionary()
sage: D.basic_variables()
(x3, x4)

c_B()

Return the $$c_B$$ vector, objective coefficients of basic variables.

OUTPUT:

• a vector

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.revised_dictionary(1, 2)
sage: D.c_B()
(10, 5)

c_N()

Return the $$c_N$$ vector, objective coefficients of non-basic variables.

OUTPUT:

• a vector

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.revised_dictionary()
sage: D.c_N()
(10, 5)

column_coefficients(v)

Return the coefficients of a nonbasic variable.

INPUT:

• v – a nonbasic variable of self, can be given as a string, an actual variable, or an integer interpreted as the index of a variable

OUTPUT:

• a vector

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.revised_dictionary()
sage: D.column_coefficients(1)
(1, 3)

constant_terms()

Return constant terms in the relations of self.

OUTPUT:

• a vector

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.revised_dictionary()
sage: D.constant_terms()
(1000, 1500)

dictionary()

Return a regular LP dictionary matching self.

OUTPUT:

EXAMPLES:

sage: A = ([1, 1], [3, 1], [-1, -1])
sage: b = (1000, 1500, -400)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.revised_dictionary()
sage: D.dictionary()
LP problem dictionary (use typeset mode to see details)

nonbasic_indices()

Return the non-basic indices of self.

Note

Non-basic indices are indices of nonbasic_variables() in the list of generators of the coordinate_ring() of the problem() of self, they may not coincide with the indices of variables which are parts of their names. (They will for the default indexed names.)

OUTPUT:

• a list

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.revised_dictionary()
sage: D.nonbasic_indices()
[1, 2]

nonbasic_variables()

Return non-basic variables of self.

OUTPUT:

• a vector

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.revised_dictionary()
sage: D.nonbasic_variables()
(x1, x2)

objective_coefficients()

Return coefficients of the objective of self.

OUTPUT:

• a vector

These are coefficients of non-basic variables when basic variables are eliminated.

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.revised_dictionary()
sage: D.objective_coefficients()
(10, 5)

objective_name()

Return the objective name of self.

OUTPUT:

• a symbolic expression

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.revised_dictionary()
sage: D.objective_name()
z

objective_value()

Return the value of the objective at the basic solution of self.

OUTPUT:

• a number

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.revised_dictionary()
sage: D.objective_value()
0

problem()

Return the original problem.

OUTPUT:

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.revised_dictionary()
sage: D.problem() is P
True

row_coefficients(v)

Return the coefficients of the basic variable v.

These are the coefficients with which nonbasic variables are subtracted in the relation for v.

INPUT:

• v – a basic variable of self, can be given as a string, an actual variable, or an integer interpreted as the index of a variable

OUTPUT:

• a vector of coefficients of a basic variable

EXAMPLES:

sage: A = ([-1, 1], [8, 2])
sage: b = (2, 17)
sage: c = (55/10, 21/10)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.revised_dictionary()
sage: D.row_coefficients("x3")
(-1, 1)

We can also use indices of variables:

sage: D.row_coefficients(3)
(-1, 1)

Or variable names without quotes after injecting them:

sage: P.inject_variables()
Defining x0, x1, x2, x3, x4
sage: D.row_coefficients(x3)
(-1, 1)

update()

Update self using previously set entering and leaving variables.

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.revised_dictionary()
sage: D.objective_value()
0
sage: D.enter("x1")
sage: D.leave("x4")
sage: D.update()
sage: D.objective_value()
5000

x_B()

Return the basic variables of self.

OUTPUT:

• a vector

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.revised_dictionary()
sage: D.basic_variables()
(x3, x4)

x_N()

Return non-basic variables of self.

OUTPUT:

• a vector

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.revised_dictionary()
sage: D.nonbasic_variables()
(x1, x2)

y()

Return the $$y$$ vector, the product of c_B() and B_inverse().

OUTPUT:

• a vector

EXAMPLES:

sage: A = ([1, 1], [3, 1])
sage: b = (1000, 1500)
sage: c = (10, 5)
sage: P = InteractiveLPProblemStandardForm(A, b, c)
sage: D = P.revised_dictionary()
sage: D.y()
(0, 0)

sage.numerical.interactive_simplex_method.default_variable_name(variable)

Return default variable name for the current style().

INPUT:

• variable - a string describing requested name

OUTPUT:

• a string with the requested name for current style

EXAMPLES:

sage: sage.numerical.interactive_simplex_method.default_variable_name("primal slack")
'x'
sage: sage.numerical.interactive_simplex_method.style('Vanderbei')
'Vanderbei'
sage: sage.numerical.interactive_simplex_method.default_variable_name("primal slack")
'w'
sage: sage.numerical.interactive_simplex_method.style('UAlberta')
'UAlberta'

sage.numerical.interactive_simplex_method.random_dictionary(m, n, bound=5, special_probability=0.2)

Construct a random dictionary.

INPUT:

• m – the number of constraints/basic variables
• n – the number of decision/non-basic variables
• bound – (default: 5) a bound on dictionary entries
• special_probability – (default: 0.2) probability of constructing a potentially infeasible or potentially optimal dictionary

OUTPUT:

EXAMPLES:

sage: from sage.numerical.interactive_simplex_method \
....:     import random_dictionary
sage: random_dictionary(3, 4)
LP problem dictionary (use typeset mode to see details)

sage.numerical.interactive_simplex_method.style(new_style=None)

Set or get the current style of problems and dictionaries.

INPUT:

• new_style – a string or None (default)

OUTPUT:

• a string with current style (same as new_style if it was given)

If the input is not recognized as a valid style, a ValueError exception is raised.

Currently supported styles are:

• ‘UAlberta’ (default): Follows the style used in the Math 373 course on Mathematical Programming and Optimization at the University of Alberta, Edmonton, Canada; based on Chvatal’s book.

• Objective functions of dictionaries are printed at the bottom.

Variable names default to

• $$z$$ for primal objective
• $$z$$ for dual objective
• $$w$$ for auxiliary objective
• $$x_1, x_2, \dots, x_n$$ for primal decision variables
• $$x_{n+1}, x_{n+2}, \dots, x_{n+m}$$ for primal slack variables
• $$y_1, y_2, \dots, y_m$$ for dual decision variables
• $$y_{m+1}, y_{m+2}, \dots, y_{m+n}$$ for dual slack variables
• ‘Vanderbei’: Follows the style of Robert Vanderbei’s textbook, Linear Programming – Foundations and Extensions.

• Objective functions of dictionaries are printed at the top.

Variable names default to

• $$zeta$$ for primal objective
• $$xi$$ for dual objective
• $$xi$$ for auxiliary objective
• $$x_1, x_2, \dots, x_n$$ for primal decision variables
• $$w_1, w_2, \dots, w_m$$ for primal slack variables
• $$y_1, y_2, \dots, y_m$$ for dual decision variables
• $$z_1, z_2, \dots, z_n$$ for dual slack variables

EXAMPLES:

sage: sage.numerical.interactive_simplex_method.style()
'UAlberta'
sage: sage.numerical.interactive_simplex_method.style('Vanderbei')
'Vanderbei'
sage: sage.numerical.interactive_simplex_method.style('Doesntexist')
Traceback (most recent call last):
...
ValueError: Style must be one of: UAlberta, Vanderbei
sage: sage.numerical.interactive_simplex_method.style('UAlberta')
'UAlberta'

sage.numerical.interactive_simplex_method.variable(R, v)

Interpret v as a variable of R.

INPUT:

• R – a polynomial ring
• v – a variable of R or convertible into R, a string with the name of a variable of R or an index of a variable in R

OUTPUT:

• a variable of R

EXAMPLES:

sage: from sage.numerical.interactive_simplex_method \
....:     import variable
sage: R = PolynomialRing(QQ, "x3, y5, x5, y")
sage: R.inject_variables()
Defining x3, y5, x5, y
sage: variable(R, "x3")
x3
sage: variable(R, x3)
x3
sage: variable(R, 3)
x3
sage: variable(R, 0)
Traceback (most recent call last):
...
ValueError: there is no variable with the given index
sage: variable(R, 5)
Traceback (most recent call last):
...
ValueError: the given index is ambiguous
sage: variable(R, 2 * x3)
Traceback (most recent call last):
...
ValueError: cannot interpret given data as a variable
sage: variable(R, "z")
Traceback (most recent call last):
...
ValueError: cannot interpret given data as a variable


text:
## Book Recommendation for Differential Forms

If anyone is seeking a more formal treatment of differential forms than the (admittedly informal!) description given in class, a good reference is

Abraham, Marsden, Ratiu, “Manifolds, Tensor Analysis, and Applications

Note that an electronic version of this book is available for free for CMU students through the library webpage.

A big difference from the treatment we’ve seen in class is that this book first spends several chapters defining and studying manifolds before introducing differential forms. We instead started with differential forms in $\mathbb{R}^n$, and will later talk about how to work with them on curves and surfaces. Interestingly enough, however, differential forms in $\mathbb{R}^n$ is essentially all we need to define discrete differential forms, which in turn are sufficiet to work with “curved” polyhedral surfaces. (The joys of being piecewise Euclidean…)

## Assignment -1: Favorite Formula

Part of your course grade is determined by participation, which can include both in-class participation as well as discussion here on the course webpage.  Therefore, your first assignment is to:

1. create an account (you must use your Andrew email address, so we can give you participation credit this semester!),
To make things interesting, your comment should include a description of your favorite mathematical formula typeset in $\LaTeX$.  If you don’t know how to use $\LaTeX$ this is a great opportunity to learn — a very basic introduction can be found here.  (And if you don’t have a favorite mathematical formula, this is a great time to pick one!)

format: plain


text:
Does perfectly normal $\implies$ normal?

A question here on perfectly normal spaces got me into investigating the definition of such a space. The definition on wikipedia says

A perfectly normal space is a topological space X in which every two disjoint non-empty closed sets $E$ and $F$ can be precisely separated by a continuous function $f$ from $X$ to the real line $\mathbb{R}$: the preimages of $\{0\}$ and $\{1\}$ under $f$ are, respectively, $E$ and $F$. (In this definition, the real line can be replaced with the unit interval $[0,1]$.) It turns out that $X$ is perfectly normal if and only if $X$ is normal and every closed set is a $G_\delta$ set.

I took it for granted that a perfectly normal space was also assumed to be a normal space, which seems to be the case if you take the definition to be a perfect space where every closed set is $G_\delta$. If you take the other definition, I don't think you can assume automatically that the space is perfect. So this is kind of a follow up, from the first definition above, how can you conclude that a perfectly normal space is also normal?

-
Related question: math.stackexchange.com/questions/72138/… –  Martin Sleziak Oct 13 '11 at 8:08
Let $F$ and $G$ be disjoint closed sets in a perfectly normal space $X$ as defined in Wikipedia. There is a continuous function $f:X\to [0,1]$ such that $F = f^{-1}[\{0\}]$ and $G = f^{-1}[\{1\}]$. Let $$U = f^{-1}\left[\left[0,\frac13\right)\right]$$ and $$V = f^{-1}\left[\left(\frac23,1\right]\right].$$ Since $f$ is continuous, $U$ and $V$ are open in $X$. Clearly $U \cap V = \varnothing$, $F\subseteq U$, and $G\subseteq V$, so $U$ and $V$ separate the closed sets $F$ and $G$. It follows that $X$ is normal.
One last thing to ask, do you know why $\mathbb{R}$ can be replaced by $[0,1]$ in this definition? –  Danielle Intal Oct 15 '11 at 4:51
@Danielle: The existence of separating functions from $X$ to $\mathbb{R}$ implies that $X$ is normal and that every closed set in $X$ is $f^{-1}[\{0\}]$ for some continuous $f:X\to\mathbb{R}$. If $F=f^{-1}[\{0\}]$, let $g(x)=\min\{|f(x)|,1\}$; $F=g^{-1}[\{0\}]$, and $g:X\to[0,1]$. Now suppose that $F,G$ are disjoint closed sets, $f,g:X\to[0,1]$ are continuous, $F=f^{-1}[\{0\}]$, and $G=g^{-1}[\{0\}]$. Let $h(x)=\dfrac{f(x)}{f(x)+g(x)}$; then $h:X\to[0,1]$ is continuous, $F=h^{-1}[\{0\}]$, and $G=h^{-1}[\{1\}]$. –  Brian M. Scott Oct 15 '11 at 9:44
I had that same $h$ in mind, but was unsure how to know $f(x)+g(x)\neq 0$ for some $x$ unless the codomain was $[0,1]$ instead of $\mathbb{R}$. Thanks for your comment! –  Danielle Intal Oct 15 '11 at 17:20


text:
# zbMATH — the first resource for mathematics

$$L^p$$ eigenfunction bounds for the Hermite operator. (English) Zbl 1075.35020
The authors obtain $$L^{p}$$ eigenfunction bounds for the Hermite operator $$H=-\Delta+x^{2}$$ in $${\mathbb R}^{n}$$ and also for a larger class of related operator of the form $$H_{V}=-\Delta+V$$. The main results are the following.
Theorem 1. (i) The operators $$e^{itH}$$ satisfy$$\| e^{H}\|_{L^{1}\to L^{\infty}}\lesssim |\sin t|^{-n/2}.$$ (ii) Let $$V$$ be a potential that satisfies $$|\partial^{\alpha}V|\leq c_{\alpha}$$, $$|\alpha|\geq 2.$$ Then for $$|t|\ll 1$$, the operators $$e^{H}_{V}$$ satisfy $$\| e^{H}_{V}\|_{L^{1}\to L^{\infty}}\lesssim |t|^{-n/2}.$$
Theorem 2. (Strichartz estimate). Let $$V$$ be a potential that satisfies $$|\partial^{\alpha}V|\lesssim 2,\quad |\alpha|=1.$$
Then the solution $$u$$ to $$(i\partial_{t}-H_{V})u=f$$, $$u(0)=u_{0},$$ satisfies $\| u\|_{L^{p_{1}}(0,1;L^{q_{1}})}\lesssim\| u_{0}\|_{L^{2}}+\| f\|_{L^{p_{2}^{\prime}}(0,1;L^{q_{2}^{\prime}})},$ whenever the pairs $$(p_{1},q_{1})$$ and $$(p_{2},q_{2})$$ are subject to $\frac{2}{p}+\frac{n-1}{q}=\frac{n-1}{2},\quad2\leq p\leq\infty,\quad(n,p,q)\neq (2,2,\infty).$ Put $$D^{\text{int}}_{j}=\{|x|\in[\lambda(1-2^{-2(j-1)}),\lambda(1-2^{-2(j+1)})]\},\quad 1\leq 2^{j}\leq\lambda^{2/3}$$, $$D^{\text{bd}}=\{||x|-\lambda|\leq\lambda^{-1/3}\},$$ $$D^{\text{ext}}=\{|x|>\lambda+\lambda^{-1/3}/2\}$$, $\| f\|_{l^{\infty}_{\lambda}L^{p}}^{q}=\| f\|^{q}_{L^{p}(D^{\text{ext}})}+\| f\|^{q}_{L^{p}(D^{\text{bd}})}+\sum_{1\leq j}^{2^{j}\leq\lambda^{2/3}}\| f\|^{q}_{L^{p}(D^{\text{int}}_{j})},$
$y=\lambda^{-2/3}(\lambda^{2}-x^{2}),\quad \langle y\rangle_{-}=1+y_{-},\quad \langle y\rangle_{+}=1+y_{+}.$
Theorem 3 (Weighted $$L^{p}$$ eigenfunction bound). (i) Let $$2\leq p\leq 2(n+1)/(n-1)$$. Then $\|\lambda^{1/3-(n/3)(1/2-1/p)}\langle y\rangle_{+}^{-1/4+((n+3)/4)(1/2-a/p)}\langle y\rangle_{-}^{1-(n/2)(1/2-1/p)}\phi\|_{l^{\infty}_{\lambda}L^{p}} \lesssim \|\phi\|_{L^{2}}+\| (H-\lambda^{2})\phi\|_{L^{2}}.$
(ii) Let $$2(n+1)/(n-1)\leq p\leq\infty$$. Let $$N$$ be a positive integer. Then $\|\lambda^{1/3-n/3(1/2-1/p)}\langle y\rangle_{+}^{1/2-(n/2)(1/2-1/p)}\langle y\rangle_{-}^{N}P_{\lambda^{2}}\phi\|_{l^{\infty}_{\lambda}L^{p}}\lesssim\|\phi\|_{L^{2}},$ where $$P_{k}$$ is a spectral projector for $$H$$.
The authors get similar estimates for $$H_{V}$$ with a positive potential $$V$$ that satisfies $$V\sim |x|^{2}$$, $$|\nabla V|\sim |x|$$, and $$|\partial^{2}_{x}V|\lesssim 1$$.

##### MSC:
 35P05 General topics in linear spectral theory for PDEs 35B60 Continuation and prolongation of solutions to PDEs 35S05 Pseudodifferential operators as generalizations of partial differential operators
Full Text:
##### References:
 [1] G. E. Andrews, R. Askey, and R. Roy, Special Functions , Encyclopedia Math. Appl. 71 , Cambridge Univ. Press, Cambridge, 1999. · Zbl 0920.33001 [2] L. Escauriaza, Carleman inequalities and the heat operator , Duke Math. J. 104 (2000), 113–127. · Zbl 0979.35029 [3] L. Escauriaza and L. Vega, Carleman inequalities and the heat operator, II , Indiana Univ. Math. J. 50 (2001), 1149–1169. · Zbl 1029.35046 [4] M. V. Fedoryuk, Asymptotic Analysis: Linear Ordinary Differential Equations , Springer, Berlin, 1993. · Zbl 0782.34001 [5] J. Ginibre and G. Velo, Generalized Strichartz inequalities for the wave equation , J. Funct. Anal. 133 (1995), 50–68. · Zbl 0849.35064 [6] D. Jerison, Carleman inequalities for the Dirac and Laplace operators and unique continuation , Adv. Math. 62 (1986), 118–134. · Zbl 0627.35008 [7] D. Jerison and C. E. Kenig, Unique continuation and absence of positive eigenvalues for Schrödinger operators , with an appendix by E. M. Stein, Ann. of Math. (2) 121 (1985), 463–494. JSTOR: · Zbl 0593.35119 [8] G. E. Karadzhov, Riesz summability of multiple Hermite series in $$L^ p$$ spaces, C. R. Acad. Bulgare Sci. 47 (1994), 5–8. · Zbl 0829.40003 [9] M. Keel and T. Tao, Endpoint Strichartz estimates , Amer. J. Math. 120 (1998), 955–980. · Zbl 0922.35028 [10] H. Koch and D. Tataru, Carleman estimates and unique continuation for second-order elliptic equations with nonsmooth coefficients , Comm. Pure Appl. Math. 54 (2001), 339–360. · Zbl 1033.35025 [11] –. –. –. –., Dispersive estimates for principally normal pseudodifferential operators , Comm. Pure Appl. Math. 58 (2005), 217–284. · Zbl 1078.35143 [12] H. Lindblad and C. D. Sogge, On existence and scattering with minimal regularity for semilinear wave equations , J. Funct. Anal. 130 (1995), 357–426. · Zbl 0846.35085 [13] H. F. Smith and C. D. Sogge, On the critical semilinear wave equation outside convex obstacles , J. Amer. Math. Soc. 8 (1995), 879–916. JSTOR: · Zbl 0860.35081 [14] C. D. Sogge, Strong uniqueness theorems for second order elliptic differential equations , Amer. J. Math. 112 (1990), 943–984. JSTOR: · Zbl 0734.35012 [15] ——–, Fourier Integrals in Classical Analysis , Cambridge Tracts in Math. 105 , Cambridge Univ. Press, Cambridge, 1993. [16] S. Thangavelu, Summability of Hermite expansions, I, II , Trans. Amer. Math. Soc. 314 (1989), 119–142., 143–170. JSTOR: · Zbl 0685.42015 [17] ——–, Lectures on Hermite and Laguerre Expansions , with a preface by R. S. Strichartz, Math. Notes 42 , Princeton Univ. Press, Princeton, 1993. · Zbl 0791.41030 [18] –. –. –. –., Hermite and special Hermite expansions revisited , Duke Math. J. 94 (1998), 257–278. · Zbl 0945.42014
This reference list is based on information provided by the publisher or from digital mathematics libraries. Its items are heuristically matched to zbMATH identifiers and may contain data conversion errors. It attempts to reflect the references listed in the original paper as accurately as possible without claiming the completeness or perfect precision of the matching.


text:
Implementing Truncated Normal Initializer

@alexis-jacq This approximation is not close to the true distribution.

Well it seems ok with plots:

in blue, the normal distribution (sorted)
in red, the truncated distribution (also sorted)

2 Likes

Well, the method confines the random variables to a circle when you want a square [Edit: The square would not work either as it has the wrong marginals, too. One needs the full “other dimension” for each coordinate]. It probably works just fine for the initialisation heuristic, but it is far from general purpose.
(Interestingly, a similar transformation to polar coordinates is the easiest way I know to compute the normalisation constant in the normal transform. If you could apply it to squares you would get a closed form for the normal CDF.)

Best regards

Thomas

I don’t understand why do you want a square… Don’t we want to confine into a circle of ray 2? And anyway, here my example is in 1d, so the circle and the square are the same, it’s the segment [-2, 2]. The box-muller method for truncated normal is already known and used, since this paper: http://oa.upm.es/22659/1/INVE_MEM_2012_153251.pdf

@alexis-jacq I haven’t looked at the paper. But from my simulation, the sampled points are not following the true truncated normal distribution,

The red line is the pdf given by scipy.stats.truncnorm

2 Likes

This seems to be the picture I had in mind when I said that it was with cutting off too much by restricting to the circle…

If you can do with an approximation, torch.fmod(torch.randn(size),2) might be easiest.
It moves more of the cut-off mass to the center than to the boundaries, but at least the good samples stay where they belong…

7 Likes

Wow, and your approximation takes one line!

So it looks like you take all the element over 2 and put them close to the mean, like in a snake game, and that’s why you have still too many people close to the mean. And it seems that the worst case comes with variance = 1. I would worth to write down the mathematics, to see how far we are from a Gaussian with both methods

Hi Alexis,

Well, I think one of the things is that your sketch seems different from the usual definition of truncated Gaussian.
So what the .fmod in sampling does in terms of pdf is to sum the pdf over all x+sign(x)*cutoff*i | i=0,1,2,… . If the shape of pdf in [cutoff*i, cutoff*(i+1)] were proportional to [0,cutoff], this would be exact. As it is falls of much more quickly than proportional (a 3 std event conditional on the magnitude being at least 2 std is much less likely than a 1 std event). This is why summation puts too much mass in the center. With a smaller cutoff, the pdf in [-cutoff,cutoff] becomes flatter (closer to uniform) and thus the error is more pronounced. And this would further increase if you cut of at less than 1 std, but at some point you might just use a uniform distribution anyway. Even at 1 you could use a cutoff of 2/3 (or so) and then multiply the sample by 3/2 to approximate the truncated normal better (i.e. torch.fmod(torch.randn(size),2/3)*3/2).

In the plain Box-Muller method, you are essentially using two transforms: first that for U_1 uniform on [0,1] R = -sqrt(2 ln(U_1)) has the distribution of the radii in a 2d standard normal, χ(2). Thus that X,Y = sin/cos(U_2) * R will be 2d standard normal. Then you are using that the marginals of the 2d standard normal are 1d standard normals.

If you now restrict R to R <= c, and take the pdf of X close to c, you have only very little of the total mass along the vertical (Y-coordinate) lies in the disc R<=c. This is why the “restricted BM” distribution has the density going to 0 as X approaches c.

But then I would be surprised if the initialisation heuristic broke down on with a simple approximation.

Best regards

Thomas

1 Like

Ok thanks, now I understand better. So, close to x= 0, all samples should be accepted, since cos(U_2) is small, but with my approach, I still multiply by a restricted R and I make it even smaller. That explains the plot of @ruotianluo where I pass over expected density close to the mean.

It was too beautiful to work…

Once again, I should be more careful when emotions mix with mathematics…

Coming back to this thread after one year.

One simple approximated implementation for truncated_normal to 2*std.

def truncated_normal_(tensor, mean=0, std=1):
size = tensor.shape
tmp = tensor.new_empty(size + (4,)).normal_()
valid = (tmp < 2) & (tmp > -2)
ind = valid.max(-1, keepdim=True)[1]
tensor.data.copy_(tmp.gather(-1, ind).squeeze(-1))

8 Likes

Hi! Does it implement same thing as https://www.tensorflow.org/api_docs/python/tf/truncated_normal ?

Here is a simpler way to sample values from truncated normal distribution.

from scipy.stats import truncnorm
import torch

def truncated_normal(size, threshold=1):
values = truncnorm.rvs(-threshold, threshold, size=size)
return values

# usage example
x= truncnorm([10, 20], threshold=1)   # sample 10x20 sized tensor
x = torch.from_numpy(x).cuda()

5 Likes

Thanks @ruotianluo.
I confirmed that your function is equivalent to scipy truncnorm.

if anyone wants reproduce:

import torch
from scipy.stats import truncnorm
import matplotlib.pyplot as plt

def truncated_normal_(tensor, mean=0, std=1):
size = tensor.shape
tmp = tensor.new_empty(size + (4,)).normal_()
valid = (tmp < 2) & (tmp > -2)
ind = valid.max(-1, keepdim=True)[1]
tensor.data.copy_(tmp.gather(-1, ind).squeeze(-1))
return tensor

fig, ax = plt.subplots(1, 1)

def test_truncnorm():
a, b = -2, 2
size = 1000000
r = truncnorm.rvs(a, b, size=size)
ax.hist(r, density=True, histtype='stepfilled', alpha=0.2, bins=50)

tensor = torch.zeros(size)
utils.truncated_normal_(tensor)
r = tensor.numpy()

ax.hist(r, density=True, histtype='stepfilled', alpha=0.2, bins=50)
ax.legend(loc='best', frameon=False)
plt.show()

if __name__ == '__main__':
test_truncnorm()

The Tensorflow implementation of truncated_normal can be found on GitHub and references https://people.sc.fsu.edu/~jburkardt/presentations/truncated_normal.pdf (see page 24) there.

Analogous Numpy/PyTorch code would look like

import numpy as np
import torch

def parameterized_truncated_normal(uniform, mu, sigma, a, b):
normal = torch.distributions.normal.Normal(0, 1)

alpha = (a - mu) / sigma
beta = (b - mu) / sigma

alpha_normal_cdf = normal.cdf(alpha)
p = alpha_normal_cdf + (normal.cdf(beta) - alpha_normal_cdf) * uniform

p = p.numpy()
one = np.array(1, dtype=p.dtype)
epsilon = np.array(np.finfo(p.dtype).eps, dtype=p.dtype)
v = np.clip(2 * p - 1, -one + epsilon, one - epsilon)
x = mu + sigma * np.sqrt(2) * torch.erfinv(torch.from_numpy(v))
x = torch.clamp(x, a, b)

return x

def truncated_normal(uniform):
return parameterized_truncated_normal(uniform, mu=0.0, sigma=1.0, a=-2, b=2)

def sample_truncated_normal(shape=()):
return truncted_normal(torch.from_numpy(np.random.uniform(0, 1, shape)))

However in the case of mu=0.0, sigma=1.0, a=-2, b=2, I suspect repeated normal sampling is also fine: 5% of entries need to be resampled each iteration. E.g. since 0.05**5 * 1000000 < 1 one would need roughly 5 iterations for 1M entries.

4 Likes

In case anyone needs a truncated normal distribution class, implementing the torch.distributions.Distribution interface, I put together one here: https://github.com/toshas/torch_truncnorm

2 Likes

Hello @anton, Thanks for posting your take on this. The class was very intuitive to use. But I found that when supplied a value outside of [a,b], log_prob does not give -inf. For example, a distribution with a=0, b=inf, loc=0, scale=1, log_prob(-1) gives -0.7. might be worth looking into.

I admit, this is not the most efficient way of doing it but it is re-drawing the values rather than clipping them as it is described in the tf specification (“values more than two standard deviations from the mean are discarded and re-drawn”)

def truncated_normal(t, mean=0.0, std=0.01):
torch.nn.init.normal_(t, mean=mean, std=std)
while True:
cond = torch.logical_or(t < mean - 2*std, t > mean + 2*std)
if not torch.sum(cond):
break
t = torch.where(cond, torch.nn.init.normal_(torch.ones(t.shape), mean=mean, std=std), t)
return t

I am using it like this m.weight.data = truncated_normal(m.weight.data) to initialize my weights. It is convenient in my case

2 Likes

Use torch.nn.init.trunc_normal_.

Description as given Here:

Fills the input Tensor with values drawn from a truncated
normal distribution. The values are effectively drawn from the
normal distribution :math:\mathcal{N}(\text{mean}, \text{std}^2)
with values outside :math:[a, b] redrawn until they are within
the bounds. The method used for generating the random values works
best when :math:a \leq \text{mean} \leq b.

1 Like

It seems torch.nn.init.trunc_normal_ does not appear in the documentation of TORCH.NN.INIT, so I am a little confused whether it is stable version of this method?

2 Likes


text:
# Netzschleuder network catalogue, repository and centrifuge

Problems with this dataset? Open an issue.
You may also take a look at the source code.
The network in this dataset can be loaded directly from graph-tool with:
import graph_tool.all as gt
g = gt.collection.ns["bookcrossing"]

# bookcrossing — BookCrossing ratings (2005)

Description

Two bipartite networks representing people and the books they have interacted with, from the BookCrossing website. Nodes represent users and books, and an edge connects a user to a book they have interacted with. The file book_implicit is unweighted; edge weights in book_ratings give the rating a user assigned to a book.1

1. Description obtained from the ICON project.

Tags
Economic Preferences Unweighted Weighted
Citation
Upstream URL OK
http://konect.cc/networks/bookcrossing_full-rating
Networks
Tip: hover your mouse over a table header to obtain a legend.
Name Nodes Edges $\left<k\right>$ $\sigma_k$ $\lambda_h$ $\tau$ $r$ $c$ $\oslash$ $S$ Kind Mode NPs EPs gt GraphML GML csv
bookcrossing 445,801 1,149,739 5.16 45.45 134.05 187.56 -0.08 0.00 19 0.94 Undirected Bipartite 9.1 MiB 13.1 MiB 12.9 MiB 11.8 MiB
Ridiculograms