Non-negative matrix factorization

In developing the hyperspectral browser project I have been playing a bit with various methods of dimension reduction, including principal components analysis and non-negative matrix factorization. Here is the base algorithm for NMF. For my application this gives interesting results with a randomly seeded input, but the solutions it finds are always local minima. It looks like there are some interesting methods for adding constraints to this algorithm, but I haven't yet been able to figure out any of the more advanced algorithms. If anyone has experience with these, please let me know! 

Here is the simple version. The input is a 3D wave, and the output components describe the layers dimension of the input data.

// the rows of matrix H contain the estimated components
// the output matrix NMFcomponents maps the proportions of components
function NMF(wave w3D, variable components)
	variable rows = DimSize(w3D, 0)
	variable cols = DimSize(w3D, 1)
	variable layers = DimSize(w3D, 2)
	variable pixels = rows * cols
	variable mincycles = 1000, maxcycles = 2000
	
	Duplicate/free w3D w_input
	Redimension/E=1/N=(pixels, layers) w_input
	
	w_input = max(0, w_input)
	
	// rows of W are pixels of w3D, columns are z-dimension of w3D
	Make/O/N=(pixels, components)/Y=(WaveType(w_input)) W
	Make/O/N=(components, layers)/Y=(WaveType(w_input)) H
	W = abs(gnoise(1))
	H = abs(gnoise(1))
	
	int i
	variable err_old = Inf
	variable start = ticks
	
	for (i=0;i<maxcycles;i++)
				
		MatrixOP/O W = W * ((w_input x H^t) / ((W x H) x H^t))
		MatrixOP/O H = H * ((W^t x w_input) / ((W^t x W) x H))
		
		if ((ticks-start)/60>10)
			Make/N=(rows, cols, components)/O NMFcomponents
			NMFcomponents = W[p + q * rows][r]
			DoUpdate
			DoAlert 1, num2str(i) + " iterations completed.\r\rContinue?"
			if (v_flag == 2)
				break
			endif
			start = ticks
		endif
				
		if (i >= mincycles)
			MatrixOP/O err = sqrt(sum(sq(w_input-(W x H))))
			if(err[0][0]/err_old > 0.9995)
				break
			endif
			err_old = err[0][0]
		endif
	endfor
	
	Make/N=(rows, cols, components)/O NMFcomponents
	NMFcomponents = W[p + q * rows][r]
	
	return i
end