Efficiently finding/ sampling from all solutions to a constrained linear problem

Question

Start with $N>3$ vectors $\vec{v}_I$ in $\mathbb{R}^3_+$, any $3$ of which are linearly independent. $I$ here ranges from $0$ to $N-1$.
Let $v_{\left[abc\right]}$ be a matrix in $\mathbb{R}^{3 \times 3}_+$ that picks three of the (column) vectors $\vec{v}_a\, \vec{v}_b\, \vec{v}_c$ in order $a < b < c$, without duplication. There are $\frac{1}{3!}N \left(N-1\right) \left(N-2\right)$ such matrices, and $\left[abc\right]$ enumerates them.

For each $v_{\left[abc\right]}$, find the linear combination $\vec{w}_{\left[abc\right]}$ that lands on a particular arbitrarily chosen point $\vec{p}$, i.e. $v_{\left[abc\right]} \vec{w}_{\left[abc\right]} = \vec{p}$.
Discard those linear combinations, which contain negative coordinates. Of the remaining, now non-negative combinations, if there are any, find the average. (Otherwise, return nothing)

$$\begin{align*} {w_{\left[abc\right]}}_i &= \sum_j {{v}_{\left[abc\right]}^{-1}}_{ij} p_j \\ \omega_{\left[abc\right]} &= \begin{cases} 1 & \text{if } 0 \le \min_i\left({w_{\left[abc\right]}}_i\right) \\ 0 & \text{else} \end{cases} \\ n &= \sum_{\left[abc\right]}\omega_{\left[abc\right]} \\ w_i &= \frac{1}{n} \sum_{\left[abc\right]} {w_{\left[abc\right]}}_i \omega_{\left[abc\right]} \text{ if } n>0 \end{align*}$$

An alternative filtering method also filters out linear combinations greater than some constant, say, $1$:

$$ \omega_{\left[abc\right]}= \begin{cases} 1 & \text{if } 0 \le \min_i\left({w_{\left[abc\right]}}_i\right) \le 1 \\ 0 & \text{else} \end{cases} $$

and instead of the unweighted average, a sample from an $n$-choice dirichlet distribution could also be taken as a weighted average, to sample from the space of all weight combinations that satisfy the filtering constraint. (I hope this is a correct statement. I haven't actually formally proved this part)

Now the question is, is there a method to calculate this, that is more efficient, than explicitly following this recipe? The $\frac{1}{3!}N \left(N-1\right) \left(N-2\right)$ very quickly makes this intractable for large $N$. In my use case, I'm looking for an $N$ somewhere between $400$ and $500$, giving me, in the worst case, $0 \le \left[abc\right] < 20708500$, i.e. over 20 million matrix inversions for evaluating the correct weights for a single point.

Note: Constrained Optimization methods can absolutely find one solution that fulfills the constraints, but this method, if I did it right, finds, and can sample from the space of all of them. This is a desirable property for me. If the ability to sample from the full space makes this difficult, just finding precisely the (unweighted) average solution in a less resource-hungry way would also be appreciated.

---

I don't need absolutely insane speeds. It just needs to be tractable. I plan to implement this in python/numpy, so I don't expect ultra-optimized peak performance. But any way to save time and memory is good. Particularly great would be a "local" version of this, that lets me directly calculate a particular $w_i$ with minimal effort.

---

I hope the notation here is clear enough. Some of the dimensionalities here are tricky so I'm not 100% sure that I got it right. Please ask if anything doesn't quite make sense.

Answer 1

I suggest you randomly sample 10,000 matrices $v_{[abc]}$. Keep only the ones where the resulting linear combination $w_{\left[abc\right]}$ meets your filtering constraints (e.g., is all positive). Average the ones kept. Output that as an estimate of the value you are seeking.

This should be a pretty good estimator of the value you are computing. You are trying to compute the expected value, with the expectation taken over all $N(N-1)(N-2)/6$ combinations, and this is a finite-sample estimate of that expectation. To put it another way, the sample mean is a good estimator for the population mean, and this method computes the sample mean. We're essentially using rejection sampling here to address the filtering condition.

How good will this be? That will depend on how many matrices/combinations you keep, with this method. Heuristically, I expect that the linear combination has about a probability $1/2$ to be non-negative in each coordinate, or about a $1/8$ probability to be non-negative in all three coordinates. So if your filtering condition is "all three coordinates must be non-negative", then I expect that you'll keep about $10000/8 \approx 1250$ combinations, which should be enough to get a good estimate of mean.

In comparison, if you had a filtering condition that was satisfied by only a miniscule fraction of matrices, then this method would not be a good choice.