Implementing NNLS - NNKL solvers#
In this lession, your goal is to
implement a NNLS/NNKL solver of your choice (AS, HALS, MU),
test is numerically on synthetic examples,
implement a NMF solver using alternating optimization.
You can start by choosing one type of loss function, either least squares
or Kullback-Leibler
Then choose one of the three algorithms seen in class:
MU is simpler, but slower for NNLS
HALS is moderately hard to implement, and fast for matrix problems
Active set is the hardest to implement (but also fun !)
Keep in mind that you will use this function with a matrix input when solving NNF later on.
Please refrain from generating all the code with AI tools (these algorithms are available online so it will copy paste the solution in seconds), but feel free to use it for debugging if you are stuck.
Packages you need: numpy (pip install numpy) tensorly (pip install tensorly)
import numpy as np
import matplotlib.pyplot as plt
def MultiplicativeUpdatesLS(y,U, x0=None,maxiter=1000):
"""
Multiplicative updates for least squares problem.
Parameters:
y : array-like, shape (n_samples,)
The observed data vector.
U : array-like, shape (n_samples, n_features)
The design matrix.
maxiter : int
Maximum number of iterations.
Returns:
x : array-like, shape (n_features,)
The estimated coefficients.
loss: array-like, shape (maxiter,)
The loss at each iteration.
"""
return x, loss
To test your algorithm, generate random nonnegative vectors \(x\) and nonnegative linear operators \(U\), then create synthetic examples of observations \(y\) according to a noise model that depends on the loss you chose.
Remember to introduce the hyperparameters explicitely (problem dimensions, noise level, number of iterations) so that you can easily play with them.
A classical metric for measuring error is the mean squared error
Also look at the loss function across iterations to debug your solver!
Implementation of a NMF solver#
You can now use your NNLS/NNKL solver to implement an algorithm for computing an approximate NMF.
Keep in mind that iterative algorithms need initialization, and this is particularly important for NMF as it is nonconvex. A classical solution is to use the absolute value of SVD outputs. But you can also use separable NMF (code provided)
Another thing consider is the number of inner iterations of the NNLS/NNKL solvers. You can set this as an hyperparameter and compare various choices on a synthetic dataset.
from auxiliary_fun import spa
def NMF(Y, r, init=None, maxiter=1000, tol=1e-6, inner_maxiter=10):
"""
Non-negative Matrix Factorization using SPA for initialization and Multiplicative Updates for optimization.
Parameters:
Y : array-like, shape (n_samples, n_features)
The data matrix to be factorized.
r : int
The rank of the factorization.
init : tuple of array-like, optional
Initial values for W and H. If None, SPA is used for initialization.
maxiter : int
Maximum number of iterations.
Returns:
W : array-like, shape (n_samples, r)
The basis matrix.
H : array-like, shape (r, n_features)
The coefficient matrix.
loss: array-like, shape (maxiter,)
The loss at each iteration.
"""
return W, H, loss
To estimate the performance of your NMF algorithm, you can generate synthetic data as in the NNLS/NNKL problem. The issue is with measuring the error for \(W\) and \(H\) which suffer from scaling and permutation ambiguitees.
You can propose a solution to fix these ambiguitees, and compare your estimated solutions with the true ones. Because NMF may not be unique, mind that these solution may not match.
On a real example#
We may now compute NMF on a real data example, anticipating on next lecture !
We will transcribe a simple piano recording, detailled explanations will be given tomorrow.
This requires installing two other packages: soundfile and scipy
import soundfile as sf
import scipy.signal as signal
from auxiliary_fun import permute_spectra
song, sr = sf.read('./MAPS_MUS-muss_1_ENSTDkAm.wav') # loads data
song = song / np.max(np.abs(song)) # normalize the audio
song = song[:, 0] # take left channel only
t_1, t_2 = 0, 6 * sr # first six seconds: isolated notes only
song_simple = song[t_1:t_2] # cut the first 6 seconds of the song for analysis
frequencies, times, Sxx = signal.stft(song_simple, fs=sr, nperseg=4096, nfft=4096, noverlap=4096 - 882) # compute stft
half_freqs = len(frequencies) // 2 # keep right side of the spectrum only
frequencies = frequencies[:half_freqs]
Sxx = Sxx[:half_freqs, :]
Sxx_abs = np.abs(Sxx)
# hyperparameters
rank = 7
maxiter = 200
W_init = np.abs(np.random.randn(Sxx.shape[0], rank))
H_init = np.abs(np.random.randn(Sxx.shape[1], rank)).T
# insert your algorithm here
W_kl, H_kl, loss = NMF(Sxx_abs, rank, init=(W_init, H_init), maxiter=100)
# Normalize W columnwise, put the norm in H, order by increasing fundamental frequency
W_kl_norms = np.max(W_kl, axis=0)
W_kl = W_kl / W_kl_norms
H_kl = H_kl.T * W_kl_norms
# Insert your postprocessing here
#perm = permute_spectra(W_kl)
#W_kl, H_kl = W_kl[:, perm], H_kl[:, perm]
# You can try and guess which note is which
#notes = ["F3", "G3", "Bb3", "Hammer", "C4", "D4", "F4"]
plt.figure(figsize=(10, 6))
for i in range(rank):
plt.subplot(rank, 2, 2 * i + 1)
plt.plot(frequencies, W_kl[:, i])
#plt.ylabel(notes[i])
if i == 0:
plt.title('Frequency profile ($W$)')
if i == rank - 1:
plt.xlabel('Frequency (Hz)')
plt.subplot(rank, 2, 2 * i + 2)
plt.plot(times, H_kl[:, i])
plt.ylim(0, np.max(H_kl))
if i == 0:
plt.title('Time activation ($H$)')
if i == rank - 1:
plt.xlabel('Time (s)')
plt.tight_layout()
plt.figure()
plt.semilogy(loss)
plt.title('Loss over iterations for NMF')
plt.xlabel('Iteration')
plt.ylabel('Loss')
plt.show()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[3], line 23
19 W_init = np.abs(np.random.randn(Sxx.shape[0], rank))
20 H_init = np.abs(np.random.randn(Sxx.shape[1], rank)).T
21
22 # insert your algorithm here
---> 23 W_kl, H_kl, loss = NMF(Sxx_abs, rank, init=(W_init, H_init), maxiter=100)
24
25 # Normalize W columnwise, put the norm in H, order by increasing fundamental frequency
26 W_kl_norms = np.max(W_kl, axis=0)
Cell In[2], line 26, in NMF(Y, r, init, maxiter, tol, inner_maxiter)
22 loss: array-like, shape (maxiter,)
23 The loss at each iteration.
24 """
25
---> 26 return W, H, loss
NameError: name 'W' is not defined
If you have finished all the coding, you can try implementing another NNLS/NNKL solver and comparing them!