Note

Go to the end to download the full example code.

Example with custom penalty class for unimodality for all but one component

In this example, we first demonstrate how to specify exactly how the penalties are imposed in the AO-ADMM fitting procedure. Then, we create a custom penalty that imposes non-negativity on all component vectors and unimodality on all but one of the component vectors.

Imports

import matplotlib.pyplot as plt
import numpy as np
import scipy.stats as stats
import tensorly as tl
from tlviz.factor_tools import factor_match_score

import matcouply.decomposition as decomposition
from matcouply.coupled_matrices import CoupledMatrixFactorization

Setup

I, J, K = 10, 40, 15
rank = 3
noise_level = 0.2
rng = np.random.default_rng(0)


def normalize(x):
    return x / tl.sqrt(tl.sum(x**2, axis=0, keepdims=True))

Generate simulated data that follows the PARAFAC2 constraint

We start by generating some components, for the $\mathbf{A}$ and $\mathbf{C}$ components, we use uniformly distributed component vector elements. For the $\mathbf{B}^{(i)}$-components, we create two unimodal vectors and one component vector with uniformly distributed elements, and shift these vectors for each $i$.

# Random uniform components
A = rng.uniform(size=(I, rank)) + 0.1  # Add 0.1 to ensure that there is signal for all components for all slices
A = tl.tensor(A)

B_0 = tl.zeros((J, rank))
# Simulating unimodal components
t = np.linspace(-10, 10, J)
for r in range(rank - 1):
    sigma = rng.uniform(0.5, 1.5)
    mu = rng.uniform(-10, 0)
    B_0[:, r] = stats.norm.pdf(t, loc=mu, scale=sigma)
# The final component is random uniform, not unimodal
B_0[:, rank - 1] = rng.uniform(size=(J,))

# Shift the components for each slice
B_is = [np.roll(B_0, i, axis=0) for i in range(I)]
B_is = [tl.tensor(B_i) for B_i in B_is]


# Random uniform components
C = rng.uniform(size=(K, rank))
C = tl.tensor(C)

Plot the simulated components

fig, axes = plt.subplots(2, 3, tight_layout=True)

axes[0, 0].plot(normalize(A))
axes[0, 0].set_title("$\\mathbf{A}$")

axes[0, 1].plot(normalize(C))
axes[0, 1].set_title("$\\mathbf{C}$")

axes[0, 2].axis("off")

axes[1, 0].plot(normalize(B_is[0]))
axes[1, 0].set_title("$\\mathbf{B}_0$")

axes[1, 1].plot(normalize(B_is[I // 2]))
axes[1, 1].set_title(f"$\\mathbf{{B}}_{{{I//2}}}$")

axes[1, 2].plot(normalize(B_is[-1]))
axes[1, 2].set_title(f"$\\mathbf{{B}}_{{{I-1}}}$")
fig.legend(["Component 0", "Component 1", "Component 2"], bbox_to_anchor=(0.95, 0.75), loc="center right")
fig.suptitle("Simulated components")

plt.show()

$Simulated components, $\mathbf{A}$, $\mathbf{C}$, $\mathbf{B}_0$, $\mathbf{B}_{5}$, $\mathbf{B}_{9}$$

For the $\mathbf{B}^{(i)}$-s, we see that component 0 and 1 are unimodal, while component 2 is not.

Create the coupled matrix factorization, simulated data matrices and add noise

cmf = CoupledMatrixFactorization((None, (A, B_is, C)))
matrices = cmf.to_matrices()
noise = [tl.tensor(rng.uniform(size=M.shape)) for M in matrices]
noisy_matrices = [M + N * noise_level * tl.norm(M) / tl.norm(N) for M, N in zip(matrices, noise)]

Use the `regs` parameter to input regularization classes

Matcouply automatically parses the constraints from the parafac2_aoadmm and cmf_aoadmm funciton arguments. However, sometimes, you may want full control over how a penalty is implemented. In that case, the regs-argument is useful. This argument makes it possible to specify exactly which penalty instances to use.

Since the components are non-negative, it makes sense to fit a non-negative PARAFAC2 model, however, we also know that two of the $\mathbf{B}^{(i)}$-component vectors are unimodal, so we first try with a fully unimodal decomposition.

from matcouply.penalties import NonNegativity, Unimodality

lowest_error = float("inf")
for init in range(4):
    print("Init:", init)
    out = decomposition.parafac2_aoadmm(
        noisy_matrices,
        rank,
        n_iter_max=1000,
        regs=[[NonNegativity()], [Unimodality(non_negativity=True)], [NonNegativity()]],
        return_errors=True,
        random_state=init,
        verbose=50,  # Only print every 50 iteration
    )
    if out[1].regularized_loss[-1] < lowest_error and out[1].satisfied_stopping_condition:
        out_cmf, diagnostics = out
        lowest_error = diagnostics.rec_errors[-1]

print("=" * 50)
print(f"Final reconstruction error: {lowest_error:.3f}")
print(f"Feasibility gap for A: {diagnostics.feasibility_gaps[-1][0]}")
print(f"Feasibility gap for B_is: {diagnostics.feasibility_gaps[-1][1]}")
print(f"Feasibility gap for C: {diagnostics.feasibility_gaps[-1][2]}")

Init: 0
All regularization penalties (including regs list):
* Mode 0:
   - <'matcouply.penalties.NonNegativity' with aux_init='random_uniform', dual_init='random_uniform')>
* Mode 1:
   - <'matcouply.penalties.Parafac2' with svd='truncated_svd', update_basis_matrices=True, update_coordinate_matrix=True, n_iter=1, aux_init='random_uniform', dual_init='random_uniform')>
   - <'matcouply.penalties.Unimodality' with non_negativity=True, aux_init='random_uniform', dual_init='random_uniform')>
* Mode 2:
   - <'matcouply.penalties.NonNegativity' with aux_init='random_uniform', dual_init='random_uniform')>
Feasibility gaps for A: [0.6536201371726535]
Feasibility gaps for the Bi-matrices: [0.991541154718251, 0.7209249811534285]
Feasibility gaps for C: [0.698917129490007]
Coupled matrix factorization iteration=0, reconstruction error=0.3314035749797261, regularized loss=0.05491416475467147 regularized loss variation=0.9067503005970164.
Feasibility gaps for A: [0.007604248629946072]
Feasibility gaps for the Bi-matrices: [0.49058444441575305, 0.216200837230665]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=50, reconstruction error=0.2988241200229145, regularized loss=0.0446479273537346 regularized loss variation=0.0041730430460376056.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.005599748831554825, 0.07667375386621993]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=100, reconstruction error=0.29785648881423815, regularized loss=0.04435924396437319 regularized loss variation=0.0011881407930202255.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0050094444203326125, 0.07180217320018303]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=150, reconstruction error=0.294138000996816, regularized loss=0.04325858181520146 regularized loss variation=0.0036106456523279746.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0034708038996178864, 0.07184365518506455]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=200, reconstruction error=0.29064079635262563, regularized loss=0.0422360362522442 regularized loss variation=0.006487353542753333.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0012741643853616687, 0.04740122551753342]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=250, reconstruction error=0.2894397507992607, regularized loss=0.04188768467136907 regularized loss variation=0.00620999116383777.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0006535342383363287, 0.02978511759271819]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=300, reconstruction error=0.2889316012131705, regularized loss=0.0417407350898033 regularized loss variation=0.0049455419494807765.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0006393959965885595, 0.023528987537592663]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=350, reconstruction error=0.28829963601881414, regularized loss=0.04155834006429036 regularized loss variation=0.005982919305416168.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.00045946963789559377, 0.016428340281644693]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=400, reconstruction error=0.2883317467091476, regularized loss=0.04156759808017402 regularized loss variation=0.0033198151656059207.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.00032441720970429247, 0.011851148272540428]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=450, reconstruction error=0.2883962258944587, regularized loss=0.04158619155508382 regularized loss variation=0.0051179503038934595.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.00034245312871417554, 0.010356397804950212]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=500, reconstruction error=0.28801236795811114, regularized loss=0.0414755620484192 regularized loss variation=0.004413430778265683.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.00019943164444400518, 0.007295906159875355]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=550, reconstruction error=0.28745181764455685, regularized loss=0.041314273733579786 regularized loss variation=0.00927009164852328.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.00023002428324823808, 0.006620289322964176]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=600, reconstruction error=0.28817402162912703, regularized loss=0.041522133370952286 regularized loss variation=0.003202841718348099.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.00015090801583741544, 0.005146935692336633]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=650, reconstruction error=0.2882077559184747, regularized loss=0.041531855285781544 regularized loss variation=0.0052045183149595515.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.00016502926783164547, 0.004844781998492417]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=700, reconstruction error=0.28767134066016153, regularized loss=0.041377400118607355 regularized loss variation=0.005534934380260884.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.00012361008089695148, 0.007050310433727602]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=750, reconstruction error=0.28833520981282507, regularized loss=0.04156859660890293 regularized loss variation=0.0006546385848659998.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.00013205904189391308, 0.004003453585765906]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=800, reconstruction error=0.28811365302016556, regularized loss=0.04150473852831218 regularized loss variation=0.0033405026221441365.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0001139054049192399, 0.0032805137953516006]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=850, reconstruction error=0.28804911449844994, regularized loss=0.04148614618167056 regularized loss variation=0.006068518277065778.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.00010613331409324228, 0.0031024189959399447]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=900, reconstruction error=0.28714333912287787, regularized loss=0.04122564860131802 regularized loss variation=0.007963031512931242.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0005503475581183409, 0.0027588118695728854]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=950, reconstruction error=0.288211569972069, regularized loss=0.04153295453288241 regularized loss variation=0.0012041085480238808.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0001731640185059626, 0.00263732164971147]
Feasibility gaps for C: [0.0]
REACHED MAXIMUM NUMBER OF ITERATIONS
Init: 1
All regularization penalties (including regs list):
* Mode 0:
   - <'matcouply.penalties.NonNegativity' with aux_init='random_uniform', dual_init='random_uniform')>
* Mode 1:
   - <'matcouply.penalties.Parafac2' with svd='truncated_svd', update_basis_matrices=True, update_coordinate_matrix=True, n_iter=1, aux_init='random_uniform', dual_init='random_uniform')>
   - <'matcouply.penalties.Unimodality' with non_negativity=True, aux_init='random_uniform', dual_init='random_uniform')>
* Mode 2:
   - <'matcouply.penalties.NonNegativity' with aux_init='random_uniform', dual_init='random_uniform')>
Feasibility gaps for A: [0.8420701798435941]
Feasibility gaps for the Bi-matrices: [0.9737546725688532, 0.691901967686556]
Feasibility gaps for C: [0.7433855015485626]
Coupled matrix factorization iteration=0, reconstruction error=0.344985210037103, regularized loss=0.05950739757217203 regularized loss variation=0.837156903159741.
Feasibility gaps for A: [0.00999730042828633]
Feasibility gaps for the Bi-matrices: [0.20663569486247232, 0.21316798877616946]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=50, reconstruction error=0.303002936596462, regularized loss=0.045905389793039784 regularized loss variation=0.0003873349555075787.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0008513823031411778, 0.0008772327964806217]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=100, reconstruction error=0.30225051736801073, regularized loss=0.045677687624615076 regularized loss variation=6.670798922240489e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [6.668298424650233e-05, 5.083663270520775e-05]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=150, reconstruction error=0.302206222335921, regularized loss=0.04566430040927405 regularized loss variation=3.0263908334107034e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.4967275268847624e-05, 1.4205180473710084e-05]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=200, reconstruction error=0.3021764997407128, regularized loss=0.0456553184977745 regularized loss variation=1.1125916087218849e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.512459639871819e-05, 3.10097544570829e-05]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=250, reconstruction error=0.3021186234043597, regularized loss=0.04563783130387266 regularized loss variation=7.136904697890996e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [5.2022630775217755e-05, 4.563359533424667e-05]
Feasibility gaps for C: [0.0]
converged in 296 iterations: FEASIBILITY GAP CRITERION AND RELATIVE LOSS CRITERION SATISFIED
Init: 2
All regularization penalties (including regs list):
* Mode 0:
   - <'matcouply.penalties.NonNegativity' with aux_init='random_uniform', dual_init='random_uniform')>
* Mode 1:
   - <'matcouply.penalties.Parafac2' with svd='truncated_svd', update_basis_matrices=True, update_coordinate_matrix=True, n_iter=1, aux_init='random_uniform', dual_init='random_uniform')>
   - <'matcouply.penalties.Unimodality' with non_negativity=True, aux_init='random_uniform', dual_init='random_uniform')>
* Mode 2:
   - <'matcouply.penalties.NonNegativity' with aux_init='random_uniform', dual_init='random_uniform')>
Feasibility gaps for A: [0.914047539381028]
Feasibility gaps for the Bi-matrices: [0.9819293155819506, 0.6953992142081842]
Feasibility gaps for C: [0.7549569841061697]
Coupled matrix factorization iteration=0, reconstruction error=0.32829205733818073, regularized loss=0.05388783745566767 regularized loss variation=0.8507497360905348.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.38228852648132033, 0.23073647629374164]
Feasibility gaps for C: [0.0009314483755839048]
Coupled matrix factorization iteration=50, reconstruction error=0.31236830260755627, regularized loss=0.04878697823696292 regularized loss variation=0.006813974457738301.
Feasibility gaps for A: [1.0624933570549639e-05]
Feasibility gaps for the Bi-matrices: [0.004992585559639334, 0.055370055664303384]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=100, reconstruction error=0.3119336065052966, regularized loss=0.04865128743370061 regularized loss variation=0.005804161385186894.
Feasibility gaps for A: [8.270516437549411e-06]
Feasibility gaps for the Bi-matrices: [0.0034403594488342927, 0.047167333037834325]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=150, reconstruction error=0.31182068055500745, regularized loss=0.048616068410894 regularized loss variation=0.005605183845065307.
Feasibility gaps for A: [6.239341621293822e-06]
Feasibility gaps for the Bi-matrices: [0.00254877144067615, 0.04327512747482581]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=200, reconstruction error=0.31166684343499906, regularized loss=0.04856811064836811 regularized loss variation=0.005913499828938771.
Feasibility gaps for A: [4.752931640083787e-06]
Feasibility gaps for the Bi-matrices: [0.001997031602371156, 0.04107598244417696]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=250, reconstruction error=0.31162049252133395, regularized loss=0.04855366567961937 regularized loss variation=0.00501859060808464.
Feasibility gaps for A: [3.7733997928188757e-06]
Feasibility gaps for the Bi-matrices: [0.0017128224561442763, 0.040541435113761945]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=300, reconstruction error=0.3113851845456219, regularized loss=0.04848036657725551 regularized loss variation=0.0038150247167217516.
Feasibility gaps for A: [3.3675121451316355e-06]
Feasibility gaps for the Bi-matrices: [0.0013984916097211637, 0.0392841929108175]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=350, reconstruction error=0.3113209041411293, regularized loss=0.04846035267762511 regularized loss variation=0.005458365600897878.
Feasibility gaps for A: [2.1920479712748834e-06]
Feasibility gaps for the Bi-matrices: [0.0012215913812582454, 0.039050846794016734]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=400, reconstruction error=0.3110550421867212, regularized loss=0.04837761963489145 regularized loss variation=0.0022709188779997495.
Feasibility gaps for A: [2.146302121446659e-06]
Feasibility gaps for the Bi-matrices: [0.0010863750845884262, 0.03833256770073923]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=450, reconstruction error=0.3111549612164023, regularized loss=0.04840870494479041 regularized loss variation=0.003127656385255317.
Feasibility gaps for A: [1.788111247941916e-06]
Feasibility gaps for the Bi-matrices: [0.0010410325127373837, 0.03826218069518964]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=500, reconstruction error=0.31091983570898324, regularized loss=0.04833557211865056 regularized loss variation=0.002581446382643314.
Feasibility gaps for A: [1.5045328206117047e-06]
Feasibility gaps for the Bi-matrices: [0.0009764126549764174, 0.037948145768054525]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=550, reconstruction error=0.31106937322688133, regularized loss=0.0483820774798824 regularized loss variation=0.003321014993838743.
Feasibility gaps for A: [1.3022301154749617e-06]
Feasibility gaps for the Bi-matrices: [0.0009567921989919036, 0.03793094673148455]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=600, reconstruction error=0.3108559579931646, regularized loss=0.04831571330992406 regularized loss variation=0.0027050062949275488.
Feasibility gaps for A: [1.1320479560733844e-06]
Feasibility gaps for the Bi-matrices: [0.0009217250949402963, 0.0377616176528814]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=650, reconstruction error=0.3110232586340214, regularized loss=0.04836773370566269 regularized loss variation=0.0034104324022380572.
Feasibility gaps for A: [1.0047060161550192e-06]
Feasibility gaps for the Bi-matrices: [0.0009125304803460744, 0.03775997158162996]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=700, reconstruction error=0.31081834077137405, regularized loss=0.048304020479935 regularized loss variation=0.002768274543791526.
Feasibility gaps for A: [8.936164345967351e-07]
Feasibility gaps for the Bi-matrices: [0.0008916682259104059, 0.037657814108479584]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=750, reconstruction error=0.31099440714584237, regularized loss=0.04835876063799698 regularized loss variation=0.003460055924219541.
Feasibility gaps for A: [8.075330410059222e-07]
Feasibility gaps for the Bi-matrices: [0.0008871280615393304, 0.0376610014395763]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=800, reconstruction error=0.3107935465899636, regularized loss=0.04829631430098394 regularized loss variation=0.002805660903330708.
Feasibility gaps for A: [7.303411821211971e-07]
Feasibility gaps for the Bi-matrices: [0.0008737596804719988, 0.03759406644322183]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=850, reconstruction error=0.31097470122390886, regularized loss=0.04835263240064969 regularized loss variation=0.0034908177567311107.
Feasibility gaps for A: [6.68945955951912e-07]
Feasibility gaps for the Bi-matrices: [0.0008714571877199039, 0.037598581545829374]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=900, reconstruction error=0.31077602180761355, regularized loss=0.048290867865283145 regularized loss variation=0.002829910193915325.
Feasibility gaps for A: [6.126679939511182e-07]
Feasibility gaps for the Bi-matrices: [0.0008623633343146045, 0.037551936098245337]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=950, reconstruction error=0.31096042399553026, regularized loss=0.048348192645739974 regularized loss variation=0.00351143977809402.
Feasibility gaps for A: [5.670366508405301e-07]
Feasibility gaps for the Bi-matrices: [0.0008612050870965656, 0.03755656163328753]
Feasibility gaps for C: [0.0]
REACHED MAXIMUM NUMBER OF ITERATIONS
Init: 3
All regularization penalties (including regs list):
* Mode 0:
   - <'matcouply.penalties.NonNegativity' with aux_init='random_uniform', dual_init='random_uniform')>
* Mode 1:
   - <'matcouply.penalties.Parafac2' with svd='truncated_svd', update_basis_matrices=True, update_coordinate_matrix=True, n_iter=1, aux_init='random_uniform', dual_init='random_uniform')>
   - <'matcouply.penalties.Unimodality' with non_negativity=True, aux_init='random_uniform', dual_init='random_uniform')>
* Mode 2:
   - <'matcouply.penalties.NonNegativity' with aux_init='random_uniform', dual_init='random_uniform')>
Feasibility gaps for A: [0.7493890212748296]
Feasibility gaps for the Bi-matrices: [0.9877127703723269, 0.6911383487531076]
Feasibility gaps for C: [0.6406461294636325]
Coupled matrix factorization iteration=0, reconstruction error=0.3373674559598866, regularized loss=0.056908400170423014 regularized loss variation=0.825279654989771.
Feasibility gaps for A: [0.008775850160562634]
Feasibility gaps for the Bi-matrices: [0.45834284129939057, 0.2812577162876005]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=50, reconstruction error=0.29395157991837595, regularized loss=0.04320376566825468 regularized loss variation=0.004044211601711358.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0016129311723680445, 0.05593480922354116]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=100, reconstruction error=0.29340842031307796, regularized loss=0.04304425055530791 regularized loss variation=0.0040839893877606046.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.001864057664228634, 0.051007395777813436]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=150, reconstruction error=0.2931404510511532, regularized loss=0.04296566202123678 regularized loss variation=0.003998429724352196.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0012148515386553314, 0.046771554489357756]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=200, reconstruction error=0.29307994734104326, regularized loss=0.04294792776671434 regularized loss variation=0.003936386052141022.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0011383533456214667, 0.0460353058465911]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=250, reconstruction error=0.29302641978195987, regularized loss=0.04293224134511668 regularized loss variation=0.003870503135952576.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0010761947515037038, 0.045071978106246385]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=300, reconstruction error=0.2929899080343091, regularized loss=0.042921543104976444 regularized loss variation=0.0038114040386124527.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0010133655245664383, 0.04387127525024263]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=350, reconstruction error=0.29296417215027387, regularized loss=0.04291400308184765 regularized loss variation=0.0037629584454842878.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0009525989176798124, 0.042450847178164754]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=400, reconstruction error=0.29294909977861183, regularized loss=0.042909587530549534 regularized loss variation=0.0037247879957863906.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0008934570527588793, 0.0408394615168342]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=450, reconstruction error=0.2923301123009648, regularized loss=0.04272844727894734 regularized loss variation=0.004798366449196394.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0014827875191875581, 0.04592567496199102]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=500, reconstruction error=0.2885811626607139, regularized loss=0.0416395437213047 regularized loss variation=0.00838816858855892.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0025613957285299848, 0.06203337425244252]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=550, reconstruction error=0.28743614027640074, regularized loss=0.04130976736849736 regularized loss variation=0.006051730404066422.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.001405216842005872, 0.041582322178372565]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=600, reconstruction error=0.2867571612466578, regularized loss=0.04111483476312085 regularized loss variation=0.007477802957495797.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0008192613861117057, 0.031227676788274872]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=650, reconstruction error=0.2869126348155401, regularized loss=0.04115943000839773 regularized loss variation=0.007220804238628738.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0005198054696551037, 0.023146561808333743]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=700, reconstruction error=0.2865629545972964, regularized loss=0.041059163473766085 regularized loss variation=0.007348268562698638.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.00034696406986300666, 0.016831294379539573]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=750, reconstruction error=0.2867990132971874, regularized loss=0.04112683701412013 regularized loss variation=0.007223567607003758.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0002627420826076849, 0.013158350081550372]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=800, reconstruction error=0.28642780296806974, regularized loss=0.041020443156557694 regularized loss variation=0.008701879878728227.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0005447118439624934, 0.012428505553802757]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=850, reconstruction error=0.28648972323614713, regularized loss=0.041038180759962094 regularized loss variation=0.007251520740616115.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.00018992189106783726, 0.00829660589623078]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=900, reconstruction error=0.28640062165060687, regularized loss=0.04101265804092703 regularized loss variation=0.007326390485964466.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.00033505012865594965, 0.008204618873378457]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=950, reconstruction error=0.2866152000984988, regularized loss=0.04107413646375125 regularized loss variation=0.007136617614333354.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0002553487604254786, 0.006968440398533865]
Feasibility gaps for C: [0.0]
REACHED MAXIMUM NUMBER OF ITERATIONS
==================================================
Final reconstruction error: 0.302
Feasibility gap for A: [0.0]
Feasibility gap for B_is: [1.8554475564054034e-05, 2.080595490589544e-05]
Feasibility gap for C: [0.0]

Compute factor match score to measure the accuracy of the recovered components

def get_stacked_CP_tensor(cmf):
    weights, factors = cmf
    A, B_is, C = factors

    stacked_cp_tensor = (weights, (A, np.concatenate(B_is, axis=0), C))
    return stacked_cp_tensor


fms, permutation = factor_match_score(
    get_stacked_CP_tensor(cmf), get_stacked_CP_tensor(out_cmf), consider_weights=False, return_permutation=True
)
print(f"Factor match score: {fms}")

Factor match score: 0.34042218544203945

Plot the recovered components

out_weights, (out_A, out_B_is, out_C) = out_cmf
out_A = out_A[:, permutation]
out_B_is = [out_B_i[:, permutation] for out_B_i in out_B_is]
out_C = out_C[:, permutation]

fig, axes = plt.subplots(2, 3, tight_layout=True)

axes[0, 0].plot(normalize(out_A))
axes[0, 0].set_title("$\\mathbf{A}$")

axes[0, 1].plot(normalize(out_C))
axes[0, 1].set_title("$\\mathbf{C}$")

axes[0, 2].axis("off")

axes[1, 0].plot(normalize(out_B_is[0]))
axes[1, 0].set_title("$\\mathbf{B}_0$")

axes[1, 1].plot(normalize(out_B_is[I // 2]))
axes[1, 1].set_title(f"$\\mathbf{{B}}_{{{I//2}}}$")

axes[1, 2].plot(normalize(out_B_is[-1]))
axes[1, 2].set_title(f"$\\mathbf{{B}}_{{{I-1}}}$")
fig.legend(["Component 0", "Component 1", "Component 2"], bbox_to_anchor=(0.95, 0.75), loc="center right")

fig.suptitle(r"Unimodality on the $\mathbf{B}^{(i)}$-components")
plt.show()

$Unimodality on the $\mathbf{B}^{(i)}$-components, $\mathbf{A}$, $\mathbf{C}$, $\mathbf{B}_0$, $\mathbf{B}_{5}$, $\mathbf{B}_{9}$$

We see that the $\mathbf{C}$-component vectors all follow the same pattern and that the the $\mathbf{A}$-component vectors all follow a similar pattern. This is not the case with the real, uncorrelated random, components. The $\mathbf{B}^{(i)}$-component vectors also follow a strange pattern with peaks jumping forwards and backwards, which we know are not the case with the real components either.

However, this strange behaviour is not too surprising, considering that there are only two uniomdal component vectors in the data. So this model that assumes all unimodal components might be too restrictive.

Create a custom penalty class for unimodality in all but one class

Now, we want to make a custom penalty that imposes unimodality on the first $R-1$ component vectors. Since unimodality are imposed column-wise, we know that this constraint is a matrix penalty (as opposed to a row-vector penalty or a multi-matrix penalty), so we import the MatrixPenalty-superclass from matcouply.penalties. We also know that unimodality is a hard constraint, so we import the HardConstraintMixin-class, which provides a penalty-method that always returns 0 and has an informative docstring.

from matcouply._doc_utils import (
    copy_ancestor_docstring,  # Helper decorator that makes it possible for ADMMPenalties to inherit a docstring
)
from matcouply._unimodal_regression import unimodal_regression  # The unimodal regression implementation
from matcouply.penalties import HardConstraintMixin, MatrixPenalty


class UnimodalAllExceptLast(HardConstraintMixin, MatrixPenalty):
    def __init__(self, non_negativity=False, aux_init="random_uniform", dual_init="random_uniform"):
        super().__init__(aux_init, dual_init)
        self.non_negativity = non_negativity

    @copy_ancestor_docstring
    def factor_matrix_update(self, factor_matrix, feasibility_penalty, aux):
        new_factor_matrix = tl.copy(factor_matrix)
        new_factor_matrix[:, :-1] = unimodal_regression(factor_matrix[:, :-1], non_negativity=self.non_negativity)
        if self.non_negativity:
            new_factor_matrix = tl.clip(new_factor_matrix, 0)
        return new_factor_matrix

Fit a non-negative PARAFAC2 model using the custom penalty class on the B mode

Now, we can fit a new model with the custom unimodality class

lowest_error = float("inf")
for init in range(4):
    print("Init:", init)
    out = decomposition.parafac2_aoadmm(
        noisy_matrices,
        rank,
        n_iter_max=1000,
        regs=[[NonNegativity()], [UnimodalAllExceptLast(non_negativity=True)], [NonNegativity()]],
        return_errors=True,
        random_state=init,
        verbose=50,  # Only print every 50 iteration
    )
    if out[1].regularized_loss[-1] < lowest_error and out[1].satisfied_stopping_condition:
        out_cmf, diagnostics = out
        lowest_error = diagnostics.rec_errors[-1]

print("=" * 50)
print(f"Final reconstruction error: {lowest_error:.3f}")
print(f"Feasibility gap for A: {diagnostics.feasibility_gaps[-1][0]}")
print(f"Feasibility gap for B_is: {diagnostics.feasibility_gaps[-1][1]}")
print(f"Feasibility gap for C: {diagnostics.feasibility_gaps[-1][2]}")

Init: 0
All regularization penalties (including regs list):
* Mode 0:
   - <'matcouply.penalties.NonNegativity' with aux_init='random_uniform', dual_init='random_uniform')>
* Mode 1:
   - <'matcouply.penalties.Parafac2' with svd='truncated_svd', update_basis_matrices=True, update_coordinate_matrix=True, n_iter=1, aux_init='random_uniform', dual_init='random_uniform')>
   - <'__main__.UnimodalAllExceptLast' with non_negativity=True, aux_init='random_uniform', dual_init='random_uniform')>
* Mode 2:
   - <'matcouply.penalties.NonNegativity' with aux_init='random_uniform', dual_init='random_uniform')>
Feasibility gaps for A: [0.6536201371726535]
Feasibility gaps for the Bi-matrices: [0.991541154718251, 0.7209249811534285]
Feasibility gaps for C: [0.698917129490007]
Coupled matrix factorization iteration=0, reconstruction error=0.23452515071154453, regularized loss=0.02750102315813634 regularized loss variation=0.9533005344936518.
Feasibility gaps for A: [0.006934761900188171]
Feasibility gaps for the Bi-matrices: [0.47777479376893245, 0.17115970603542544]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=50, reconstruction error=0.10750045365430454, regularized loss=0.005778173767940639 regularized loss variation=0.0013405716899638279.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.00016623857103926254, 7.246566317699339e-05]
Feasibility gaps for C: [2.6506291951745657e-05]
Coupled matrix factorization iteration=100, reconstruction error=0.10419576717356194, regularized loss=0.0054283789484435635 regularized loss variation=0.0009630728407163166.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.00016750078160197954, 3.6166358943861015e-05]
Feasibility gaps for C: [7.268217201774797e-06]
Coupled matrix factorization iteration=150, reconstruction error=0.10267973576902277, regularized loss=0.005271564068798167 regularized loss variation=0.0003645893170675969.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [6.921858448259306e-05, 1.210577533526796e-05]
Feasibility gaps for C: [1.575964802001983e-06]
Coupled matrix factorization iteration=200, reconstruction error=0.10228923375381758, regularized loss=0.005231543670971567 regularized loss variation=4.5911546507836023e-05.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.9530079003893777e-05, 3.491440592080208e-06]
Feasibility gaps for C: [4.266161848255958e-07]
Coupled matrix factorization iteration=250, reconstruction error=0.10223459686380995, regularized loss=0.00522595639795287 regularized loss variation=9.476502771305978e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [7.029351231314765e-06, 1.4496807015875652e-06]
Feasibility gaps for C: [1.4913567748125203e-07]
Coupled matrix factorization iteration=300, reconstruction error=0.10221835570407531, regularized loss=0.005224296121422433 regularized loss variation=4.7623035196974375e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [3.307500280061131e-06, 9.540021706484545e-07]
Feasibility gaps for C: [7.222596020730206e-08]
Coupled matrix factorization iteration=350, reconstruction error=0.1022111274466267, regularized loss=0.0052235572869552835 regularized loss variation=7.102901890897967e-07.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [2.671670896073069e-06, 2.090986826453026e-06]
Feasibility gaps for C: [3.492870240142188e-08]
Coupled matrix factorization iteration=400, reconstruction error=0.10220912796373646, regularized loss=0.005223352919553727 regularized loss variation=4.943398315949946e-07.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [8.440672895627757e-07, 2.5092130548915427e-07]
Feasibility gaps for C: [1.8819619870995417e-08]
Coupled matrix factorization iteration=450, reconstruction error=0.10220822705576495, regularized loss=0.005223260838941401 regularized loss variation=2.572464318371244e-07.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [4.554301975171541e-07, 1.5299990766612888e-07]
Feasibility gaps for C: [1.1117932793583882e-08]
Coupled matrix factorization iteration=500, reconstruction error=0.10220771274812471, regularized loss=0.005223208272601587 regularized loss variation=1.5943058892225237e-07.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [2.671475308703376e-07, 9.723478433073137e-08]
Feasibility gaps for C: [7.2405583808525266e-09]
Coupled matrix factorization iteration=550, reconstruction error=0.10220737567652999, regularized loss=0.005223173821341667 regularized loss variation=1.1046498677083325e-07.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.6193507051851867e-07, 6.257492869059405e-08]
Feasibility gaps for C: [4.887272823064566e-09]
Coupled matrix factorization iteration=600, reconstruction error=0.10220713227080787, regularized loss=0.005223148943511208 regularized loss variation=8.299666353460313e-08.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.0046617912434704e-07, 4.1036833803620333e-08]
Feasibility gaps for C: [3.358215778780599e-09]
Coupled matrix factorization iteration=650, reconstruction error=0.10220694385990943, regularized loss=0.005223129686591339 regularized loss variation=6.604998611039427e-08.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [6.375153340039454e-08, 2.7991341295588785e-08]
Feasibility gaps for C: [2.3210818374120335e-09]
Coupled matrix factorization iteration=700, reconstruction error=0.10220679073659217, regularized loss=0.005223114036336771 regularized loss variation=5.47187125472843e-08.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [4.182358664976404e-08, 2.064148299741619e-08]
Feasibility gaps for C: [1.5987344759600875e-09]
Coupled matrix factorization iteration=750, reconstruction error=0.1022066620032762, regularized loss=0.0052231008789259716 regularized loss variation=4.6617763353640313e-08.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [2.932120904449408e-08, 1.7034041707268076e-08]
Feasibility gaps for C: [1.0895583060619166e-09]
Coupled matrix factorization iteration=800, reconstruction error=0.10220655119384683, regularized loss=0.005223089553470217 regularized loss variation=4.0499065451189445e-08.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [2.3061778169017286e-08, 1.5569054791508153e-08]
Feasibility gaps for C: [7.32178137942494e-10]
Coupled matrix factorization iteration=850, reconstruction error=0.10220645423101668, regularized loss=0.005223079643238454 regularized loss variation=3.5667355499077796e-08.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [2.060625009140354e-08, 1.5073463065329567e-08]
Feasibility gaps for C: [4.892420527361048e-10]
Coupled matrix factorization iteration=900, reconstruction error=0.10220636840214371, regularized loss=0.00522307087097736 regularized loss variation=3.1714550636814264e-08.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.9985280336491266e-08, 1.4894580871028375e-08]
Feasibility gaps for C: [3.395183083338221e-10]
Coupled matrix factorization iteration=950, reconstruction error=0.10220629181217035, regularized loss=0.00522306304299726 regularized loss variation=2.839066900816425e-08.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.998794330136179e-08, 1.4749036100362089e-08]
Feasibility gaps for C: [2.6896861171071897e-10]
REACHED MAXIMUM NUMBER OF ITERATIONS
Init: 1
All regularization penalties (including regs list):
* Mode 0:
   - <'matcouply.penalties.NonNegativity' with aux_init='random_uniform', dual_init='random_uniform')>
* Mode 1:
   - <'matcouply.penalties.Parafac2' with svd='truncated_svd', update_basis_matrices=True, update_coordinate_matrix=True, n_iter=1, aux_init='random_uniform', dual_init='random_uniform')>
   - <'__main__.UnimodalAllExceptLast' with non_negativity=True, aux_init='random_uniform', dual_init='random_uniform')>
* Mode 2:
   - <'matcouply.penalties.NonNegativity' with aux_init='random_uniform', dual_init='random_uniform')>
Feasibility gaps for A: [0.8420701798435941]
Feasibility gaps for the Bi-matrices: [0.9737546725688532, 0.691901967686556]
Feasibility gaps for C: [0.7433855015485626]
Coupled matrix factorization iteration=0, reconstruction error=0.28228485865454495, regularized loss=0.03984237071280821 regularized loss variation=0.8909706137886058.
Feasibility gaps for A: [0.009954018187787994]
Feasibility gaps for the Bi-matrices: [0.22399754513699222, 0.17952115245757755]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=50, reconstruction error=0.0875435486741435, regularized loss=0.0038319364572310666 regularized loss variation=0.0024619847369702647.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.00032008523853734174, 0.00029937909984316136]
Feasibility gaps for C: [5.47648826181471e-06]
Coupled matrix factorization iteration=100, reconstruction error=0.08610152443153316, regularized loss=0.00370673625471695 regularized loss variation=0.00019467728802031407.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [2.3242625553125425e-05, 7.568335624250606e-06]
Feasibility gaps for C: [6.018491861031892e-07]
Coupled matrix factorization iteration=150, reconstruction error=0.08566923630679904, regularized loss=0.0036696090246950874 regularized loss variation=6.046423198739785e-05.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.3912219842641695e-05, 7.73918557338162e-06]
Feasibility gaps for C: [7.603533564159626e-07]
Coupled matrix factorization iteration=200, reconstruction error=0.08558989432012344, regularized loss=0.0036628150048649494 regularized loss variation=2.208073212961566e-05.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.7034164111595187e-05, 1.3725987132632539e-05]
Feasibility gaps for C: [4.585327868469141e-07]
Coupled matrix factorization iteration=250, reconstruction error=0.08556226678342614, regularized loss=0.003660450748559094 regularized loss variation=8.993080869557672e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [4.8331227446180925e-06, 1.67076183310371e-06]
Feasibility gaps for C: [2.5290159855859294e-07]
Coupled matrix factorization iteration=300, reconstruction error=0.08554938728821672, regularized loss=0.0036593488326946482 regularized loss variation=4.39129190943217e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [3.4858066418695386e-06, 9.81636891400477e-07]
Feasibility gaps for C: [1.687581471416167e-07]
Coupled matrix factorization iteration=350, reconstruction error=0.08554220901718726, regularized loss=0.0036587347617700766 regularized loss variation=2.574947240673477e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [2.4006338968312143e-06, 7.051699990366287e-07]
Feasibility gaps for C: [1.1048268569373337e-07]
Coupled matrix factorization iteration=400, reconstruction error=0.08553776232221365, regularized loss=0.0036583543915457563 regularized loss variation=1.747272857815927e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.6895476515421473e-06, 4.593776559840573e-07]
Feasibility gaps for C: [7.422848410617028e-08]
Coupled matrix factorization iteration=450, reconstruction error=0.08553458711919791, regularized loss=0.0036580827968258284 regularized loss variation=1.3227947468978455e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.1921430891599502e-06, 3.1996225252702774e-07]
Feasibility gaps for C: [4.8575265512299805e-08]
Coupled matrix factorization iteration=500, reconstruction error=0.08553197246357755, regularized loss=0.003657859156755094 regularized loss variation=1.1491760811277274e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [8.714455208309536e-07, 2.3382418689318717e-07]
Feasibility gaps for C: [3.391130207304488e-08]
Coupled matrix factorization iteration=550, reconstruction error=0.08553016663316909, regularized loss=0.0036577047021488355 regularized loss variation=4.6705358318471964e-07.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [6.801621562194997e-07, 3.1456735382186503e-07]
Feasibility gaps for C: [2.708840993831746e-08]
Coupled matrix factorization iteration=600, reconstruction error=0.08552942886193328, regularized loss=0.0036576416007242526 regularized loss variation=2.5468346560120575e-07.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [5.39020421101438e-07, 2.2498319472150596e-07]
Feasibility gaps for C: [2.439864846183738e-08]
Coupled matrix factorization iteration=650, reconstruction error=0.08552901056564477, regularized loss=0.003657605824169087 regularized loss variation=1.4986417728279286e-07.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [3.953926092765208e-07, 1.6535261820829499e-07]
Feasibility gaps for C: [1.9334167139324188e-08]
Coupled matrix factorization iteration=700, reconstruction error=0.08552875744488554, regularized loss=0.003657584175033032 regularized loss variation=9.36708052417248e-08.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [2.868526258868969e-07, 1.2268586364217112e-07]
Feasibility gaps for C: [1.4624891019625881e-08]
Coupled matrix factorization iteration=750, reconstruction error=0.08552859399001952, regularized loss=0.003657570194954802 regularized loss variation=6.28841271741944e-08.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [2.0708048308766122e-07, 9.12917679707028e-08]
Feasibility gaps for C: [1.078529049047223e-08]
Coupled matrix factorization iteration=800, reconstruction error=0.08552847936312306, regularized loss=0.0036575603910840838 regularized loss variation=4.6613396746582505e-08.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.4500872912897998e-07, 6.506575460048862e-08]
Feasibility gaps for C: [7.723668100263656e-09]
Coupled matrix factorization iteration=850, reconstruction error=0.0855283888283115, regularized loss=0.0036575526477834195 regularized loss variation=3.924808169527213e-08.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.0071497010253631e-07, 4.686010926496132e-08]
Feasibility gaps for C: [5.348101439175401e-09]
Coupled matrix factorization iteration=900, reconstruction error=0.08552830860229409, regularized loss=0.0036575457861846264 regularized loss variation=3.636161999580944e-08.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [7.172104614286609e-08, 3.511909054819947e-08]
Feasibility gaps for C: [3.6561134268979543e-09]
Coupled matrix factorization iteration=950, reconstruction error=0.08552823186189584, regularized loss=0.0036575392227111074 regularized loss variation=3.5672364210132675e-08.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [5.365912285421809e-08, 2.828551017212386e-08]
Feasibility gaps for C: [2.4911416469576742e-09]
REACHED MAXIMUM NUMBER OF ITERATIONS
Init: 2
All regularization penalties (including regs list):
* Mode 0:
   - <'matcouply.penalties.NonNegativity' with aux_init='random_uniform', dual_init='random_uniform')>
* Mode 1:
   - <'matcouply.penalties.Parafac2' with svd='truncated_svd', update_basis_matrices=True, update_coordinate_matrix=True, n_iter=1, aux_init='random_uniform', dual_init='random_uniform')>
   - <'__main__.UnimodalAllExceptLast' with non_negativity=True, aux_init='random_uniform', dual_init='random_uniform')>
* Mode 2:
   - <'matcouply.penalties.NonNegativity' with aux_init='random_uniform', dual_init='random_uniform')>
Feasibility gaps for A: [0.914047539381028]
Feasibility gaps for the Bi-matrices: [0.9819293155819506, 0.6953992142081842]
Feasibility gaps for C: [0.7549569841061697]
Coupled matrix factorization iteration=0, reconstruction error=0.22477174868817396, regularized loss=0.025261169504169814 regularized loss variation=0.9300354886525041.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.38539255102332187, 0.13910021179771864]
Feasibility gaps for C: [0.0004440734510898217]
Coupled matrix factorization iteration=50, reconstruction error=0.103147796470465, regularized loss=0.005319733958356236 regularized loss variation=0.016393928497415895.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.005516295343159779, 0.016255368366385282]
Feasibility gaps for C: [6.7192039536810825e-06]
Coupled matrix factorization iteration=100, reconstruction error=0.08843959760502175, regularized loss=0.003910781212269085 regularized loss variation=0.00014158123803643128.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [4.021911960776802e-05, 1.5229150308926729e-05]
Feasibility gaps for C: [2.1058715483724747e-07]
Coupled matrix factorization iteration=150, reconstruction error=0.08824914414838446, regularized loss=0.0038939557214611696 regularized loss variation=7.17406896661002e-05.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.529877575039469e-05, 6.1364419675509506e-06]
Feasibility gaps for C: [4.8563022202792456e-08]
Coupled matrix factorization iteration=200, reconstruction error=0.08816767468728692, regularized loss=0.0038867694298816277 regularized loss variation=9.307394510860008e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [3.120681126256353e-06, 1.2561851492817913e-06]
Feasibility gaps for C: [7.232388050214074e-09]
Coupled matrix factorization iteration=250, reconstruction error=0.08815442354116777, regularized loss=0.0038856011949377975 regularized loss variation=4.303969930615332e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.5976449209378455e-06, 7.452299889429851e-07]
Feasibility gaps for C: [1.5649069544889303e-09]
Coupled matrix factorization iteration=300, reconstruction error=0.08814257038057964, regularized loss=0.003884556356647718 regularized loss variation=2.1065866156289384e-05.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [4.524371082211133e-06, 1.7681147127721951e-06]
Feasibility gaps for C: [2.633490127564818e-08]
Coupled matrix factorization iteration=350, reconstruction error=0.08810919945411852, regularized loss=0.0038816155142228194 regularized loss variation=1.1350540070219445e-05.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [3.7347852606195035e-06, 1.1364499448566804e-06]
Feasibility gaps for C: [2.1508126276986266e-08]
Coupled matrix factorization iteration=400, reconstruction error=0.08808925929963296, regularized loss=0.003879858801978986 regularized loss variation=7.570921137110646e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.4464836332923594e-06, 4.624272096318e-07]
Feasibility gaps for C: [5.281461287186992e-09]
Coupled matrix factorization iteration=450, reconstruction error=0.08807334282953487, regularized loss=0.003878456858584391 regularized loss variation=7.193449440525898e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [8.123825550026755e-07, 3.5684527729928606e-07]
Feasibility gaps for C: [3.937306719945412e-09]
Coupled matrix factorization iteration=500, reconstruction error=0.0880567208921773, regularized loss=0.0038769930471414067 regularized loss variation=7.936801591588177e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [9.483883330759219e-07, 3.738483898845786e-07]
Feasibility gaps for C: [6.454721488232754e-09]
Coupled matrix factorization iteration=550, reconstruction error=0.08803801817061777, regularized loss=0.003875346321705012 regularized loss variation=9.0895733206681e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [9.945415177091302e-07, 3.4860163987414886e-07]
Feasibility gaps for C: [7.30987932119365e-09]
Coupled matrix factorization iteration=600, reconstruction error=0.08801646782983037, regularized loss=0.0038734493046197823 regularized loss variation=1.0493273751703901e-05.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.0750839306276506e-06, 3.4635104032485073e-07]
Feasibility gaps for C: [7.656780677257715e-09]
Coupled matrix factorization iteration=650, reconstruction error=0.08799181974037958, regularized loss=0.0038712801706117265 regularized loss variation=1.1887851645177452e-05.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.1447032994324334e-06, 3.535997101750877e-07]
Feasibility gaps for C: [7.75889834553505e-09]
Coupled matrix factorization iteration=700, reconstruction error=0.0879647920856407, regularized loss=0.0038689023233349986 regularized loss variation=9.049174026116803e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.4523115151522402e-06, 7.434176208232658e-07]
Feasibility gaps for C: [7.346668132237407e-09]
Coupled matrix factorization iteration=750, reconstruction error=0.08794459286038797, regularized loss=0.0038671257066897015 regularized loss variation=9.00132024831405e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.119832971937604e-06, 3.2336679102771335e-07]
Feasibility gaps for C: [5.005284891608334e-09]
Coupled matrix factorization iteration=800, reconstruction error=0.08793516958658897, regularized loss=0.0038662970251110806 regularized loss variation=8.162128526092121e-07.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [4.6615129368063037e-07, 1.3128703878204155e-07]
Feasibility gaps for C: [4.057845088477358e-09]
Coupled matrix factorization iteration=850, reconstruction error=0.08793354452160092, regularized loss=0.0038661541260661854 regularized loss variation=6.82434093734566e-07.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [2.114626529817422e-07, 7.474506917377733e-08]
Feasibility gaps for C: [1.945615419950588e-09]
Coupled matrix factorization iteration=900, reconstruction error=0.08793215525298706, regularized loss=0.00386603196371771 regularized loss variation=5.859876600137897e-07.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.4762792156182452e-07, 5.245213316499098e-08]
Feasibility gaps for C: [1.2171402919242122e-09]
Coupled matrix factorization iteration=950, reconstruction error=0.08793096213595408, regularized loss=0.003865927051077295 regularized loss variation=5.032384143820895e-07.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.0308682383669387e-07, 3.84783815571313e-08]
Feasibility gaps for C: [8.220238175499609e-10]
REACHED MAXIMUM NUMBER OF ITERATIONS
Init: 3
All regularization penalties (including regs list):
* Mode 0:
   - <'matcouply.penalties.NonNegativity' with aux_init='random_uniform', dual_init='random_uniform')>
* Mode 1:
   - <'matcouply.penalties.Parafac2' with svd='truncated_svd', update_basis_matrices=True, update_coordinate_matrix=True, n_iter=1, aux_init='random_uniform', dual_init='random_uniform')>
   - <'__main__.UnimodalAllExceptLast' with non_negativity=True, aux_init='random_uniform', dual_init='random_uniform')>
* Mode 2:
   - <'matcouply.penalties.NonNegativity' with aux_init='random_uniform', dual_init='random_uniform')>
Feasibility gaps for A: [0.7493890212748296]
Feasibility gaps for the Bi-matrices: [0.9877127703723269, 0.6911383487531076]
Feasibility gaps for C: [0.6406461294636325]
Coupled matrix factorization iteration=0, reconstruction error=0.23457041055428263, regularized loss=0.027511638753802352 regularized loss variation=0.9155336821898677.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.43631616733567324, 0.2218860628548828]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=50, reconstruction error=0.08574361603418637, regularized loss=0.003675983845308991 regularized loss variation=0.0021871892088718135.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.00020508776126804763, 4.71712060702908e-05]
Feasibility gaps for C: [8.779796213399544e-06]
Coupled matrix factorization iteration=100, reconstruction error=0.08406783885619859, regularized loss=0.0035337007649758868 regularized loss variation=0.0002469036628240984.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [5.743008602343044e-05, 9.41647259343093e-06]
Feasibility gaps for C: [1.2022907087255262e-06]
Coupled matrix factorization iteration=150, reconstruction error=0.08379989981858157, regularized loss=0.0035112116048021543 regularized loss variation=5.943296835576698e-05.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [2.125871647063141e-05, 3.545798620523273e-06]
Feasibility gaps for C: [4.108290291827059e-07]
Coupled matrix factorization iteration=200, reconstruction error=0.08372954054052642, regularized loss=0.0035053179795638285 regularized loss variation=1.7131386754467742e-05.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.0434863952556623e-05, 1.973513726504353e-06]
Feasibility gaps for C: [1.9204419517294577e-07]
Coupled matrix factorization iteration=250, reconstruction error=0.08370913304441004, regularized loss=0.00350360947752337 regularized loss variation=4.84877059056885e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [5.17685257038753e-06, 9.981002990421024e-07]
Feasibility gaps for C: [9.150305945471877e-08]
Coupled matrix factorization iteration=300, reconstruction error=0.08370336249794501, regularized loss=0.0035031264467311933 regularized loss variation=1.4307725091124693e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [2.7375564510794724e-06, 5.173019093435659e-07]
Feasibility gaps for C: [4.634201020317028e-08]
Coupled matrix factorization iteration=350, reconstruction error=0.08370168485312446, regularized loss=0.003502986023625882 regularized loss variation=3.9202217584587647e-07.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.3173624988065822e-06, 2.5232090377263923e-07]
Feasibility gaps for C: [2.2113205024703627e-08]
Coupled matrix factorization iteration=400, reconstruction error=0.08370121255617519, regularized loss=0.0035029464916870095 regularized loss variation=1.229242960125899e-07.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [6.510098400001608e-07, 1.179880804001862e-07]
Feasibility gaps for C: [1.0738626923702572e-08]
Coupled matrix factorization iteration=450, reconstruction error=0.08370105022731171, regularized loss=0.0035029329045774787 regularized loss variation=4.832060331203759e-08.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [3.2659455554410175e-07, 5.8691736724809715e-08]
Feasibility gaps for C: [5.372873837693824e-09]
Coupled matrix factorization iteration=500, reconstruction error=0.08370097936218332, regularized loss=0.003502926973094319 regularized loss variation=2.390479834701568e-08.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.6422765775509952e-07, 2.97901605139126e-08]
Feasibility gaps for C: [2.715305561247519e-09]
Coupled matrix factorization iteration=550, reconstruction error=0.08370094131155849, regularized loss=0.003502923788220479 regularized loss variation=1.39296864370174e-08.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [8.346179217592188e-08, 1.564239618197708e-08]
Feasibility gaps for C: [1.392041703385234e-09]
converged in 586 iterations: FEASIBILITY GAP CRITERION AND RELATIVE LOSS CRITERION SATISFIED
==================================================
Final reconstruction error: 0.084
Feasibility gap for A: [0.0]
Feasibility gap for B_is: [5.1439866448041737e-08, 1.0143057504688387e-08]
Feasibility gap for C: [8.65877892631776e-10]

Compute factor match score to measure the accuracy of the recovered components

fms, permutation = factor_match_score(
    get_stacked_CP_tensor(cmf), get_stacked_CP_tensor(out_cmf), consider_weights=False, return_permutation=True
)
print(f"Factor match score: {fms}")

Factor match score: 0.9750290134582585

We see that the factor match score is much better now compared to before!

Plot the recovered components

out_weights, (out_A, out_B_is, out_C) = out_cmf
out_A = out_A[:, permutation]
out_B_is = [out_B_i[:, permutation] for out_B_i in out_B_is]
out_C = out_C[:, permutation]

fig, axes = plt.subplots(2, 3, tight_layout=True)

axes[0, 0].plot(normalize(out_A))
axes[0, 0].set_title("$\\mathbf{A}$")

axes[0, 1].plot(normalize(out_C))
axes[0, 1].set_title("$\\mathbf{C}$")

axes[0, 2].axis("off")

axes[1, 0].plot(normalize(out_B_is[0]))
axes[1, 0].set_title("$\\mathbf{B}_0$")

axes[1, 1].plot(normalize(out_B_is[I // 2]))
axes[1, 1].set_title(f"$\\mathbf{{B}}_{{{I//2}}}$")

axes[1, 2].plot(normalize(out_B_is[-1]))
axes[1, 2].set_title(f"$\\mathbf{{B}}_{{{I-1}}}$")
fig.legend(["Component 0", "Component 1", "Component 2"], bbox_to_anchor=(0.95, 0.75), loc="center right")
fig.suptitle(r"Custom uniomdality on the $\mathbf{B}^{(i)}$-components")
plt.show()

$Custom uniomdality on the $\mathbf{B}^{(i)}$-components, $\mathbf{A}$, $\mathbf{C}$, $\mathbf{B}_0$, $\mathbf{B}_{5}$, $\mathbf{B}_{9}$$

We see that the model finds much more sensible component vectors. The $\mathbf{A}$- and $\mathbf{C}$-component vectors no longer seem correlated, and the peaks of the $\mathbf{B}^{(i)}$-component vectors no longer jump around.

Automated unit test creation with MatCoupLy

Matcouply can autogenerate unit tests for new penalty types. Simply create a test class that inherits from matcouply.testing.BaseTestFactorMatrixPenalty (or matcouply.testing.BaseTestFactorMatricesPenalty or matcouply.testing.BaseTestRowPenalty if your penalty type is either a row-wise penalty or a multi matrix penalty), overload the get_invariant_matrix and get_non_invariant_matrix methods (without changing the signature) and MatCoupLy will generate a set of useful unit tests. The get_invariant_matrix method should generate a matrix that is not modified by the proximal operator and the get_non_invariant_matrix should generate a matrix that will be modified by the proximal operator.

Note that to use this functionality in your own project, you need to add the MatCoupLy test fixtures to your test suite. You can do this by adding the line pytest_plugins = ["matcouply.testing.fixtures"] to your conftest.py file.

from matcouply.testing import MixinTestHardConstraint, BaseTestFactorMatrixPenalty

class TestUnimodalAllExceptLast(
    MixinTestHardConstraint, BaseTestFactorMatrixPenalty
):
    PenaltyType = UnimodalAllExceptLast
    penalty_default_kwargs = {}
    min_rows = 3
    min_columns = 2

    def get_invariant_matrix(self, rng, shape):
        matrix = tl.zeros(shape)
        I, J = shape
        t = np.linspace(-10, 10, I)
        for j in range(J-1):
            sigma = rng.uniform(0.5, 1)
            mu = rng.uniform(-5, 5)
            matrix[:, j] = stats.norm.pdf(t, loc=mu, scale=sigma)
        matrix[:, J-1] = rng.uniform(size=I)
        return matrix

    def get_non_invariant_matrix(self, rng, shape):
        # There are at least 3 rows
        M = rng.uniform(size=shape)
        M[1, :-1] = -1  # M is positive, so setting the second element to -1 makes it impossible for it to be unimodal
        return M

print("""
Auto generated unit tests:
==========================\
""")
for name in dir(TestUnimodalAllExceptLast):
    if "test" in name:
        print(f" * {name}")

Auto generated unit tests:
==========================
 * test_aux_as_matrix
 * test_auxes_as_matrices
 * test_factor_matrices_update_changes_input
 * test_factor_matrices_update_invariant_point
 * test_factor_matrices_update_reduces_penalty
 * test_factor_matrix_update_changes_input
 * test_factor_matrix_update_invariant_point
 * test_factor_matrix_update_reduces_penalty
 * test_given_init_aux
 * test_given_init_dual
 * test_input_validation_for_init_aux
 * test_input_validation_init_dual
 * test_penalty
 * test_rank_and_mode_validation_for_init_aux
 * test_rank_and_mode_validation_for_init_dual
 * test_standard_normal_init_aux
 * test_standard_normal_init_dual
 * test_subtract_from_aux
 * test_subtract_from_auxes
 * test_uniform_init_aux
 * test_uniform_init_dual
 * test_validating_given_init_aux
 * test_validating_given_init_dual
 * test_zeros_init_aux
 * test_zeros_init_dual

Total running time of the script: (1 minutes 17.260 seconds)

Gallery generated by Sphinx-Gallery

Example with custom penalty class for unimodality for all but one component

Imports

Setup

Generate simulated data that follows the PARAFAC2 constraint

Plot the simulated components

Create the coupled matrix factorization, simulated data matrices and add noise

Use the regs parameter to input regularization classes

Compute factor match score to measure the accuracy of the recovered components

Plot the recovered components

Create a custom penalty class for unimodality in all but one class

Fit a non-negative PARAFAC2 model using the custom penalty class on the B mode

Compute factor match score to measure the accuracy of the recovered components

Plot the recovered components

Automated unit test creation with MatCoupLy

Use the `regs` parameter to input regularization classes