# python version of nnls/fnnls

Seems the last email didn’t deliver the attachement (even though I received it back?)

The attachment is here again but just in case, here’s the code as plain text:

#!/usr/bin/env python

fnnls.py (4.38 KB)

···

# gjok - 20050816

import sys, math

import numpy, numpy.linalg.linalg as la

def any(X) : return len(numpy.nonzero(X)) != 0

def find(X) : return numpy.nonzero(X)

def norm(X, d) : return max(numpy.sum(abs(X)))

# x, w = fnnls(XtX, Xty, tol)

def fnnls(XtX, Xty, tol = 0) :

#FNNLS Non-negative least-squares.

# initialize variables

``````m = XtX.shape[0]

n = XtX.shape[1]

if tol == 0 :

eps = 2.2204e-16

tol = 10 * eps * norm(XtX,1)*max(m, n);

#end

P = numpy.zeros(n, numpy.Int16)

P[:] = -1

Z = numpy.arange(0,n)

z = numpy.zeros(m, numpy.Float)

x = numpy.array(P)

ZZ = numpy.array(Z)

w = Xty - numpy.dot(XtX, x)
``````

# set up iteration criterion

``````iter = 0

itmax = 30 * n
``````

# outer loop to put variables into set to hold positive coefficients

``````while any(Z) and any(w[ZZ] > tol) :

wt = w[ZZ].max()

t = find(w[ZZ] == wt)

t = t[-1:][0]

t = ZZ[t]

P[t] = t

Z[t] = -1

PP = find(P != -1)

ZZ = find(Z != -1)

if len(PP) == 1 :

XtyPP = Xty[PP]

XtXPP = XtX[PP, PP]

z[PP] = XtyPP / XtXPP

else :

XtyPP = numpy.array(Xty[PP])

XtXPP = numpy.array(XtX[PP, numpy.array(PP)[:, numpy.NewAxis]])

z[PP] = numpy.dot(XtyPP, la.generalized_inverse(XtXPP))

#end

z[ZZ] = 0
``````

# inner loop to remove elements from the positive set which no longer belong

``````    while any(z[PP] <= tol) and (iter < itmax) :

iter += 1

iztol = find(z <= tol)

ip = find(P[iztol] != -1)

QQ = iztol[ip]

if len(QQ) == 1 : alpha = x[QQ] / (x[QQ] - z[QQ])

else :

x_xz = x[QQ] / (x[QQ] - z[QQ])

alpha = x_xz.min()

#end

x += alpha * (z - x)

iabs = find(abs(x) < tol)

ip = find(P[iabs] != -1)

ij = iabs[ip]

Z[ij] = numpy.array(ij)

P[ij] = -1

PP = find(P != -1)

ZZ = find(Z != -1)

if len(PP) == 1 :

XtyPP = Xty[PP]

XtXPP = XtX[PP, PP]

z[PP] = XtyPP / XtXPP

else :

XtyPP = numpy.array(Xty[PP])

XtXPP = numpy.array(XtX[PP, numpy.array(PP)[:, numpy.NewAxis]])

z[PP] = numpy.dot(XtyPP, la.generalized_inverse(XtXPP))

#endif

z[ZZ] = 0

#end while

x = numpy.array(z)

w = Xty - numpy.dot(XtX, x)

#end while

return x, w
``````

#end def

if name == ‘main’ :

# x = lsqnonneg(X, y) => x = [0.9312; 0.3683; 0; 0];

``````X = numpy.array([[1, 10, 4, 10], [4, 5, 1, 12], [5, 1, 9, 20]], numpy.Float)

y = numpy.array([4, 7, 4], numpy.Float)

if False :

X = numpy.array([[1, 10, 4, 10],

[4, 5, 1, 12],

[5, 1, 9, 20],

[4, 3, 2, 1]], numpy.Float)

y = numpy.array([4, 7, 4, 1], numpy.Float)

#end

if False :

X = numpy.zeros((20, 20), numpy.Float)

for n in range(20) : X[n,:] = numpy.arange(0.0, 400.0, step = 20)

y = numpy.arange(0.0, 20.0)

#end

Xt = numpy.transpose(numpy.array(X))

x, w = fnnls(numpy.dot(Xt, X), numpy.dot(Xt, y))

print 'X = ', X

print 'y = ', y

print 'x = ', x
``````

#end if name == ‘main

Begin forwarded message:

From: Graeme O’Keefe <gjok@…484…>

Date: 25 May 2006 8:46:07 PM

To: matplotlib user list matplotlib-users@lists.sourceforge.net

Subject: [Matplotlib-users] Fwd: python version of nnls/fnnls

Dear matplotlibers,

I seem to use Non Negative Least Squares every couple of months and the most recent problem I solved reminded me that I meant to post this last November.

fnnls.py is a port of fnnls.m which was written by Rasmus Bro (http://newton.foodsci.kvl.dk/users/rasmus.html, code available at http://www.ub.es/gesq/mcr/als/fnnls.m, and also available on Matlab Central).

I’ve retained Rasmus’ comments and translated the matlab to python/numpy.