How to find a convex decompositions of a given element of a polytope?

Question

Suppose we are given a finite number of points $x_k\in\mathbb R^n$, and we know that $x$ is in the convex hull of these points.

Is there a procedure or algorithm to find a set of coefficients $p_k\ge0$ such that $x=\sum_k p_k x_k$ and $\sum_k p_k=1$?

I'm aware of a related question on mathoverflow asking about how to solve the associated membership problem. The answers there suggest that a possible solution comes from linear programming, but I don't know if this is still the case if we are looking for the coefficients, not just to solve membership. Also, as I'm not well-versed with linear programming, I was hoping for some more detail on how that approach would actually work (if linear programming is indeed what's needed to find the coefficients).

Answer 1

Yes, you can use linear programming to find such $p_k$. For each $k$, let $p_k$ be the desired weight. The objective is arbitrary, so you can minimize the constant $0$ function. The (linear) constraints are: \begin{align} \sum_k x^k_i p_k &= x_i &&\text{for $i \in \{1,\dots,n\}$} \\ \sum_k p_k &= 1 \\ 0 \le p_k &\le 1 &&\text{for all $k$}\\ \end{align}

Answer 2

As per the other answer, linear programming achieves this. The idea is that the problem of finding the coefficients $p_i\ge0$ such that $\sum_{i=1}^n p_i \mathbf x_i = \mathbf x\in\mathbb R^d$ with $\sum_i p_i=1$ can be reframed as the task of finding the vector $\mathbf p\in\mathbb R^n$ such that $X\mathbf p = \mathbf x,$ where $X$ is the matrix whose $i$-th column is $\mathbf x_i$, under the constraints $\mathbf p\ge0$ and $\sum_i p_i = 1$. This is a linear program in standard form.

Approaching it as a standard linear algebraic task, we can define the quantities $$ \tilde X=\begin{pmatrix} X \\ 1 \cdots 1\end{pmatrix} \quad\text{ and }\quad \tilde{\mathbf x}\equiv \begin{pmatrix}\mathbf x \\ 1\end{pmatrix},\tag1 $$ to rewrite the problem as $\tilde X\mathbf p =\tilde{\mathbf x}$. Here, $\tilde X$ is a matrix of dimensions $(d+1)\times n$. This equation admits solutions iff $\tilde{\mathbf x}\in\operatorname{range}(\tilde X)$, in which case the set of solutions has the form $$\mathbf p \in \tilde X^+\tilde{\mathbf x} + \ker(\tilde X),\tag2$$ with $\tilde X^+$ denoting the pseudo-inverse of $\tilde X$. For (2) to be an actual solution of our convex problem, we also need $p_i\ge 0$ for all $i$. We then use the shifts allowed by $\ker(\tilde X)$ to achieve this.

---

As a concrete example of this method in action, consider in $\mathbb R^2$ the points $$ \mathbf x_1 = (1, 1), \qquad \mathbf x_2 = (-1, 1), \qquad \mathbf x_3 = (-1, -1), \qquad \mathbf x_4 = (1, -1),$$ that is, the vertices of a square around the origin. Suppose we then have $\mathbf x=(0.5, 0.8)$, and wish to find the convex decomposition of $\mathbf x$ in terms of $\mathbf x_i$. We start by computing $\tilde X$ and $\tilde X^+$, which read $$ \tilde X = \begin{pmatrix} 1 & -1 & -1 & 1 \\ 1 & 1 & -1 & -1 \\ 1 & 1 & 1 & 1 \end{pmatrix}, \qquad \tilde X^+ = \frac14\begin{pmatrix} 1 & 1 & 1 \\ -1 & 1 & 1 \\ -1 & -1 & 1 \\ 1 & -1 & 1 \end{pmatrix}. $$ Applying $\tilde X^+$ to $\tilde{\mathbf x}$ we get $\tilde X^+ \tilde{\mathbf x}=(0.575, 0.325, -0.075, 0.175)^T$. One element of this vector is negative, which means that this is not a valid convex combination. We then consider the one-dimensional space obtained translating this point with the kernel of $\tilde X$. This is one-dimensional in this case: $\ker(\tilde X)=\mathbb R (-1, 1, -1, 1)^T$. The full set of solutions is therefore $$\begin{pmatrix} 0.575 \\ 0.325 \\ -0.075 \\ 0.175\end{pmatrix} + t \begin{pmatrix} -1 \\ 1 \\ -1 \\ 1\end{pmatrix}, \qquad t\in\mathbb R.$$ We can easily find that all the elements are positive for $t\in[-0.175, -0.075]$, which gives us the set of viable convex decompositions of $\mathbf x$.

Here's a visualisation of the set of solutions found with this technique:

$\hspace{100pt}$

Of course, the complexity of this naive approach might scale badly with the number of points and dimension, so for the general case optimised linear programming algorithms will perform much better.

---

A slightly more complex variation of the above is obtained asking to decompose $\mathbf x$ using five different vertices. In this case the kernel of $\tilde X$ is two-dimensional, and thus there is a two-dimensional convex region of possible coefficients. We can visualise this situation as follows:

$\hspace{100pt}$

where the figure on the left is used to choose a point in $\tilde X^+\tilde{\mathbf x}+\ker(\tilde X)$ such that all elements are positive, and on the right we show the corresponding convex decomposition of the red point in terms of the black ones.

---

The above animation was generated with Mathematica using the following code

decoupleOneElement[list_List,index_Integer]:={
    Normal @ SparseArray[index -> 1, Length @ list],
    ReplacePart[list, index -> 0] // # / Total @ # &
};

lineRepresentationConvexComb[coefficients_]:=With[{
        nonzeroIndices = Flatten @ Position[
            coefficients, _?(# !=0 &)
        ]
    },
    Fold[
        Last @ Sow @ decoupleOneElement[#1, #2]&,
        coefficients[[nonzeroIndices]],
        Range[Length @ nonzeroIndices - 1]
    ] // Reap // Last // Last // Map[
        ReplacePart[
            ConstantArray[0, Length @ coefficients],
            Thread @ Rule[nonzeroIndices, #]
        ]&,
        #, {2}
    ]&
];

DynamicModule[{pts, x, XTilde, kernelOfXTilde},
    pts = {{1, 1}, {-1, 1}, {-1, -1}, {1, -1}};
    x = {0.5, 0.8};
    XTilde = Append[Transpose@pts, ConstantArray[1, Length@pts]] // Echo[#, "X=", MatrixForm]&;
    kernelOfXTilde = First @ NullSpace @ XTilde // Echo[#,"ker(X)=", MatrixForm]&;
    Manipulate[
        Dot[PseudoInverse @ XTilde, Append[x,1]] + t * kernelOfXTilde // lineRepresentationConvexComb //
        Map[Total[# * pts]&, #, {2}] & // Map @ Line // Deploy @ Graphics[{
                PointSize @ 0.03, Point @ pts,
                {Red, Point @ x},
                #
            },
            Axes -> True, AxesOrigin -> {0, 0},
            GridLines -> Automatic,
            PlotRange -> ConstantArray[{-1.2, 1.2}, 2]
        ]&,
        {{t, -0.1}, -0.175, -0.075, 0.001, Appearance->"Labeled"}
    ]
]

---

The code to generate the second visualisation is the following:

DynamicModule[{
        pts, x, XTilde, kernelOfXTilde, solutionSet, solutionConstraints,
        pointInConstraints = {-0.1, 0.1}, movingPoint = False
    },
    pts = {{1, 1}, {0, 2}, {-1, 1}, {-1, -1}, {1, -1}};
    x = {0.5, 1};
    XTilde = Append[Transpose @ pts, ConstantArray[1, Length @ pts]] // Echo[#, "X=", MatrixForm]&;
    kernelOfXTilde = NullSpace @ XTilde // Echo[#,"ker(X)=", MatrixForm] &;
    solutionSet = Plus[
        Dot[PseudoInverse @ XTilde, Append[x,1]],
        {t1, t2} * kernelOfXTilde // Apply @ Sequence
    ];
    solutionConstraints = And @@ Thread @ GreaterEqual[solutionSet, 0];
    Row @ {
        (* plot constraints on the parameters, and allow to choose point in the feasible region *)
        EventHandler[
            RegionPlot[
                solutionConstraints,
                {t1, -4, 4}, {t2, -4, 4}, PlotRange -> All, PlotPoints -> 400,
                ImageSize -> 300, FrameLabel->{"t1", "t2"}, Frame -> True, PlotRangePadding -> None
            ] ~ Show ~ Graphics[{PointSize @ 0.03, Point @ Dynamic @ pointInConstraints}],
            {
                "MouseDown" :> If[
                    Not @ TrueQ @ movingPoint && Norm[pointInConstraints - MousePosition["Graphics"]] <= 0.1,
                    movingPoint = True
                ],
                "MouseUp" :> If[TrueQ @ movingPoint, movingPoint = False],
                "MouseDragged" :> With[{mp = MousePosition["Graphics"]},
                    If[TrueQ @ movingPoint && TrueQ[solutionConstraints /. Thread @ Rule[{t1, t2}, mp]],
                        pointInConstraints = mp
                    ]
                ]
            }
        ],
        (* show convex combination corresponding to the chosen values of t1 and t2 *)
        Dynamic @ Deploy @ Graphics[{
                PointSize @ 0.03, Point @ pts,
                {Red, Point @ x},
                solutionSet /. Thread @ Rule[{t1, t2}, pointInConstraints] //
                lineRepresentationConvexComb //
                Map[Total[# * pts]&, #, {2}] & // Map @ Line
            },
            Axes -> True, AxesOrigin -> {0, 0},
            GridLines -> Automatic, ImageSize -> 300,
            PlotRange -> All
        ]
    }
]