# Scipy.optimize.nnls usage problem in pymc3

**URL:** <https://discourse.pymc.io/t/scipy-optimize-nnls-usage-problem-in-pymc3/220>\
**Category:** Questions\
**Tags:** from\_github, theano\
**Created:** [August 4, 2017, 2:55pm UTC](https://discourse.pymc.io/t/scipy-optimize-nnls-usage-problem-in-pymc3/220 "2017-08-04T14:55:27Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![sezenyg](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/sezenyg/32/1013_2.png) [@sezenyg](https://discourse.pymc.io/u/sezenyg)\
**Post date:** [August 4, 2017, 2:55pm UTC](https://discourse.pymc.io/t/scipy-optimize-nnls-usage-problem-in-pymc3/220/1 "2017-08-04T14:55:27Z")

</div>

Hello everyone,

I am trying to implement a graphical modelling in pymc3 and I probabilistically generate a matrix and a vector. I need to find a third vector which provides A\*x=b equality where x has nonnegative elements. So, I have tried to use scipy.optimize.nnls() function but due to the pm.Model structure, my A and x matrices have zero shapes. The given sample script below provides ‘ValueError: expected matrix’ error since ‘mat’ matrix is empty. Does anyone have a solution for this problem?

```python
with Normal_model:
    mu = pm.Uniform('mu', lower=-10, upper=10, shape=5)
    sigma = pm.Uniform('sigma', lower=0, upper=10, shape=2)
    mat = pm.Normal('matrix', mu=mu[0], sd=sigma[0], shape=(5,2))
    t = pm.Deterministic('parametre', optimize.nnls(mat, mu))
    for i in range(0,2):
        data = pm.Normal('data_%i'%i, mu=mu[i], sd=sigma[i], observed=s[i, :])
    start = pm.find_MAP()
    step = pm.Metropolis()
    trace = pm.sample(10000, step, start=start, progressbar=True)

```

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [August 4, 2017, 3:02pm UTC](https://discourse.pymc.io/t/scipy-optimize-nnls-usage-problem-in-pymc3/220/2 "2017-08-04T15:02:07Z")

</div>

Did you try to warp it as a theano op using `as_op`? There are similar discussion here you might found useful: [Fitting a distribution with custom functions](https://discourse.pymc.io/t/fitting-a-distribution-with-custom-functions/210/3)

Also, just a small tip: you can wrap your code as a Markdown code broke [https://github.com/adam-p/markdown-here/wiki/Markdown-Cheatsheet#code](https://github.com/adam-p/markdown-here/wiki/Markdown-Cheatsheet#code)

---

<div class="post-metadata">

**Author:** ![aseyboldt](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/aseyboldt/32/5795_2.png) [@aseyboldt](https://discourse.pymc.io/u/aseyboldt)\
**Post date:** [August 4, 2017, 3:05pm UTC](https://discourse.pymc.io/t/scipy-optimize-nnls-usage-problem-in-pymc3/220/3 "2017-08-04T15:05:35Z")

</div>

You can also try to use `tt.nlinalg.MatrixPinv`:

```auto
t = pm.Deterministic('parametre', tt.nlinalg.MatrixPinv()(mat).dot(mu))

```

I’m guessing this is less stable than the scipy function, but it might be good enough.

---

<div class="post-metadata">

**Author:** ![sezenyg](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/sezenyg/32/1013_2.png) [@sezenyg](https://discourse.pymc.io/u/sezenyg)\
**Post date:** [August 4, 2017, 7:41pm UTC](https://discourse.pymc.io/t/scipy-optimize-nnls-usage-problem-in-pymc3/220/4 "2017-08-04T19:41:36Z")

</div>

Hi aseyboldt,

Since I need a non-negative result, I can’t only take the pseudo inverse and nnls function does what exactly I need.  
Anyways, thanks for your consideration 🙂

---

<div class="post-metadata">

**Author:** ![aseyboldt](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/aseyboldt/32/5795_2.png) [@aseyboldt](https://discourse.pymc.io/u/aseyboldt)\
**Post date:** [August 4, 2017, 8:02pm UTC](https://discourse.pymc.io/t/scipy-optimize-nnls-usage-problem-in-pymc3/220/5 "2017-08-04T20:02:07Z")

</div>

Ah, I didn’t read the docstring of `nnls` properly, and missed the `x>=0` part. If you can get by without gradient-based methods this should work:

```auto
@theano.as_op([tt.dmatrix, tt.dvector], [tt.dvector])
def nnls(mat, mu):
    return optimize.nnls(mat, mu)[0]

with pm.Model():
    mu = pm.Uniform('mu', lower=-10, upper=10, shape=5)
    sigma = pm.Uniform('sigma', lower=0, upper=10, shape=2)
    mat = pm.Normal('matrix', mu=mu[0], sd=sigma[0], shape=(5,2))
    t = pm.Deterministic('parametre', nnls(mat, mu))
    #...
    trace = pm.sample(step=pm.Metropolis())

```

If you need gradients, you’ll have to do a bit more work. [This](http://users.cecs.anu.edu.au/~sgould/papers/argmin-TR-2016.pdf) might be of help, it contains a section about inequality constraints (haven’t read that part though).

---

<div class="post-metadata">

**Author:** ![sezenyg](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/sezenyg/32/1013_2.png) [@sezenyg](https://discourse.pymc.io/u/sezenyg)\
**Post date:** [August 4, 2017, 8:30pm UTC](https://discourse.pymc.io/t/scipy-optimize-nnls-usage-problem-in-pymc3/220/6 "2017-08-04T20:30:38Z")

</div>

It works like a charm thanks a lot! 👍🏽