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.33140357497972656, regularized loss=0.05491416475467162 regularized loss variation=0.9067503005970162.
Feasibility gaps for A: [0.007604248629946068]
Feasibility gaps for the Bi-matrices: [0.49058444441575283, 0.21620083723066505]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=50, reconstruction error=0.2988241200229145, regularized loss=0.0446479273537346 regularized loss variation=0.004173043046035756.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.005599748831554883, 0.07667375386621979]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=100, reconstruction error=0.29785648881423743, regularized loss=0.04435924396437298 regularized loss variation=0.0011881407930250689.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.005009444420332623, 0.07180217320018271]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=150, reconstruction error=0.294138000996816, regularized loss=0.04325858181520146 regularized loss variation=0.0036106456523227193.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0034708038996178726, 0.0718436551850641]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=200, reconstruction error=0.2906407963526246, regularized loss=0.0422360362522439 regularized loss variation=0.006487353542760515.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0012741643853616828, 0.04740122551753343]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=250, reconstruction error=0.2894397507992607, regularized loss=0.04188768467136907 regularized loss variation=0.0062099911638358075.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0006535342383363047, 0.02978511759271829]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=300, reconstruction error=0.2889316012131705, regularized loss=0.0417407350898033 regularized loss variation=0.0049455419494845625.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0006393959965884551, 0.023528987537592524]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=350, reconstruction error=0.28829963601881536, regularized loss=0.04155834006429071 regularized loss variation=0.005982919305409188.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.00045946963789568837, 0.016428340281644485]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=400, reconstruction error=0.28833174670914885, regularized loss=0.041567598080174385 regularized loss variation=0.003319815165604067.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.00032441720970431454, 0.0118511482725399]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=450, reconstruction error=0.2883962258944584, regularized loss=0.04158619155508374 regularized loss variation=0.005117950303886698.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0003424531287141952, 0.010356397804949919]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=500, reconstruction error=0.2880123679581127, regularized loss=0.04147556204841965 regularized loss variation=0.004413430778242586.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.00019943164444399884, 0.007295906159875131]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=550, reconstruction error=0.2874518176445548, regularized loss=0.04131427373357919 regularized loss variation=0.009270091648536106.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.00023002428324844752, 0.006620289322963987]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=600, reconstruction error=0.2881740216291288, regularized loss=0.0415221333709528 regularized loss variation=0.0032028417183748007.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.00015090801583748052, 0.005146935692336504]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=650, reconstruction error=0.28820775591847575, regularized loss=0.04153185528578185 regularized loss variation=0.005204518314926115.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.00016502926783205698, 0.004844781998492402]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=700, reconstruction error=0.2876713406601559, regularized loss=0.04137740011860574 regularized loss variation=0.00553493438027901.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.00012361008089704792, 0.007050310433727181]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=750, reconstruction error=0.2883352098128147, regularized loss=0.04156859660889993 regularized loss variation=0.0006546385849977462.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.00013205904189325112, 0.004003453585765592]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=800, reconstruction error=0.28811365302017217, regularized loss=0.04150473852831408 regularized loss variation=0.0033405026220946618.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.00011390540491911186, 0.0032805137953512025]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=850, reconstruction error=0.28804911449844794, regularized loss=0.041486146181669986 regularized loss variation=0.006068518277109318.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.00010613331409357487, 0.003102418995939432]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=900, reconstruction error=0.28714333912289086, regularized loss=0.04122564860132175 regularized loss variation=0.007963031512945766.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0005503475581184049, 0.0027588118695725202]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=950, reconstruction error=0.2882115699720558, regularized loss=0.04153295453287861 regularized loss variation=0.0012041085479921713.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0001731640185059485, 0.0026373216497112104]
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.3449852100371028, regularized loss=0.05950739757217197 regularized loss variation=0.8371569031597411.
Feasibility gaps for A: [0.009997300428286312]
Feasibility gaps for the Bi-matrices: [0.2066356948624725, 0.21316798877616952]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=50, reconstruction error=0.30300293659646177, regularized loss=0.045905389793039715 regularized loss variation=0.00038733495550908963.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0008513823031413051, 0.0008772327964805842]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=100, reconstruction error=0.30225051736801045, regularized loss=0.04567768762461499 regularized loss variation=6.670798920721421e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [6.66829842463452e-05, 5.083663270515268e-05]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=150, reconstruction error=0.30220622233592126, regularized loss=0.04566430040927413 regularized loss variation=3.0263908282442777e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.4967275268891493e-05, 1.4205180473709155e-05]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=200, reconstruction error=0.30217649974071265, regularized loss=0.045655318497774454 regularized loss variation=1.1125916108496621e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.5124596398667568e-05, 3.100975445708357e-05]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=250, reconstruction error=0.30211862340435897, regularized loss=0.04563783130387244 regularized loss variation=7.136904699867556e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [5.202263077513785e-05, 4.563359533429539e-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.3822885264813205, 0.23073647629374167]
Feasibility gaps for C: [0.0009314483755838847]
Coupled matrix factorization iteration=50, reconstruction error=0.3123683026075561, regularized loss=0.04878697823696287 regularized loss variation=0.006813974457733539.
Feasibility gaps for A: [1.0624933570549691e-05]
Feasibility gaps for the Bi-matrices: [0.004992585559639296, 0.05537005566430376]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=100, reconstruction error=0.3119336065052969, regularized loss=0.0486512874337007 regularized loss variation=0.00580416138518646.
Feasibility gaps for A: [8.270516437524903e-06]
Feasibility gaps for the Bi-matrices: [0.0034403594488343577, 0.04716733303783444]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=150, reconstruction error=0.3118206805550072, regularized loss=0.04861606841089393 regularized loss variation=0.005605183845070961.
Feasibility gaps for A: [6.23934162128624e-06]
Feasibility gaps for the Bi-matrices: [0.0025487714406762277, 0.04327512747482584]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=200, reconstruction error=0.31166684343499906, regularized loss=0.04856811064836811 regularized loss variation=0.005913499828941594.
Feasibility gaps for A: [4.752931640076788e-06]
Feasibility gaps for the Bi-matrices: [0.0019970316023713267, 0.0410759824441769]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=250, reconstruction error=0.31162049252133417, regularized loss=0.04855366567961944 regularized loss variation=0.005018590608086472.
Feasibility gaps for A: [3.7733997928124086e-06]
Feasibility gaps for the Bi-matrices: [0.0017128224561441213, 0.040541435113761924]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=300, reconstruction error=0.3113851845456227, regularized loss=0.048480366577255754 regularized loss variation=0.003815024716719602.
Feasibility gaps for A: [3.3675121451347475e-06]
Feasibility gaps for the Bi-matrices: [0.001398491609720905, 0.03928419291081706]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=350, reconstruction error=0.3113209041411298, regularized loss=0.04846035267762526 regularized loss variation=0.005458365600896162.
Feasibility gaps for A: [2.192047971269368e-06]
Feasibility gaps for the Bi-matrices: [0.0012215913812583892, 0.0390508467940163]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=400, reconstruction error=0.31105504218672075, regularized loss=0.04837761963489132 regularized loss variation=0.0022709188780038964.
Feasibility gaps for A: [2.14630212144931e-06]
Feasibility gaps for the Bi-matrices: [0.001086375084587961, 0.03833256770073877]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=450, reconstruction error=0.31115496121640207, regularized loss=0.04840870494479034 regularized loss variation=0.003127656385249624.
Feasibility gaps for A: [1.7881112479419528e-06]
Feasibility gaps for the Bi-matrices: [0.0010410325127375392, 0.038262180695189135]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=500, reconstruction error=0.3109198357089837, regularized loss=0.0483355721186507 regularized loss variation=0.0025814463826325956.
Feasibility gaps for A: [1.5045328206117344e-06]
Feasibility gaps for the Bi-matrices: [0.000976412654976059, 0.037948145768054366]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=550, reconstruction error=0.3110693732268823, regularized loss=0.04838207747988269 regularized loss variation=0.0033210149938355886.
Feasibility gaps for A: [1.3022301154706793e-06]
Feasibility gaps for the Bi-matrices: [0.0009567921989920141, 0.037930946731484304]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=600, reconstruction error=0.3108559579931663, regularized loss=0.04831571330992458 regularized loss variation=0.0027050062949289483.
Feasibility gaps for A: [1.1320479560754336e-06]
Feasibility gaps for the Bi-matrices: [0.0009217250949403445, 0.03776161765288077]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=650, reconstruction error=0.3110232586340202, regularized loss=0.04836773370566231 regularized loss variation=0.0034104324022352344.
Feasibility gaps for A: [1.0047060161540504e-06]
Feasibility gaps for the Bi-matrices: [0.0009125304803450705, 0.0377599715816292]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=700, reconstruction error=0.31081834077137055, regularized loss=0.048304020479933915 regularized loss variation=0.0027682745438153016.
Feasibility gaps for A: [8.936164345902097e-07]
Feasibility gaps for the Bi-matrices: [0.000891668225910104, 0.03765781410847876]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=750, reconstruction error=0.3109944071458449, regularized loss=0.04835876063799778 regularized loss variation=0.003460055924183575.
Feasibility gaps for A: [8.075330410005462e-07]
Feasibility gaps for the Bi-matrices: [0.0008871280615396909, 0.037661001439575334]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=800, reconstruction error=0.3107935465899624, regularized loss=0.04829631430098356 regularized loss variation=0.0028056609033161573.
Feasibility gaps for A: [7.303411821194685e-07]
Feasibility gaps for the Bi-matrices: [0.0008737596804730486, 0.03759406644322067]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=850, reconstruction error=0.3109747012239098, regularized loss=0.048352632400649985 regularized loss variation=0.0034908177567156986.
Feasibility gaps for A: [6.689459559544015e-07]
Feasibility gaps for the Bi-matrices: [0.0008714571877193849, 0.037598581545827965]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=900, reconstruction error=0.31077602180761643, regularized loss=0.04829086786528404 regularized loss variation=0.0028299101939317037.
Feasibility gaps for A: [6.12667993951101e-07]
Feasibility gaps for the Bi-matrices: [0.0008623633343146464, 0.03755193609824374]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=950, reconstruction error=0.31096042399552737, regularized loss=0.04834819264573908 regularized loss variation=0.003511439778084109.
Feasibility gaps for A: [5.670366508405175e-07]
Feasibility gaps for the Bi-matrices: [0.0008612050870965273, 0.03755656163328585]
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.33736745595988615, regularized loss=0.05690840017042286 regularized loss variation=0.8252796549897715.
Feasibility gaps for A: [0.008775850160562673]
Feasibility gaps for the Bi-matrices: [0.4583428412993905, 0.2812577162876005]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=50, reconstruction error=0.2939515799183765, regularized loss=0.043203765668254845 regularized loss variation=0.004044211601715228.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.001612931172368029, 0.055934809223541375]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=100, reconstruction error=0.29340842031307796, regularized loss=0.04304425055530791 regularized loss variation=0.004083989387759142.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.001864057664228611, 0.05100739577781399]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=150, reconstruction error=0.2931404510511532, regularized loss=0.04296566202123678 regularized loss variation=0.0039984297243595215.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0012148515386551557, 0.04677155448935797]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=200, reconstruction error=0.29307994734104326, regularized loss=0.04294792776671434 regularized loss variation=0.0039363860521424866.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0011383533456215686, 0.04603530584659107]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=250, reconstruction error=0.29302641978195987, regularized loss=0.04293224134511668 regularized loss variation=0.0038705031359507844.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0010761947515036353, 0.04507197810624628]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=300, reconstruction error=0.2929899080343088, regularized loss=0.04292154310497636 regularized loss variation=0.0038114040386202793.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0010133655245664263, 0.0438712752502427]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=350, reconstruction error=0.29296417215027437, regularized loss=0.0429140030818478 regularized loss variation=0.003762958445485904.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0009525989176798709, 0.04245084717816476]
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.0008934570527588743, 0.040839461516834145]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=450, reconstruction error=0.292330112300965, regularized loss=0.0427284472789474 regularized loss variation=0.004798366449201306.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0014827875191875824, 0.04592567496199077]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=500, reconstruction error=0.28858116266071365, regularized loss=0.041639543721304635 regularized loss variation=0.00838816858855944.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0025613957285300112, 0.06203337425244218]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=550, reconstruction error=0.28743614027640046, regularized loss=0.041309767368497285 regularized loss variation=0.006051730404064563.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0014052168420058422, 0.04158232217837249]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=600, reconstruction error=0.28675716124665757, regularized loss=0.041114834763120785 regularized loss variation=0.007477802957499748.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.000819261386111717, 0.031227676788274743]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=650, reconstruction error=0.28691263481554063, regularized loss=0.04115943000839789 regularized loss variation=0.007220804238632644.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0005198054696551161, 0.02314656180833338]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=700, reconstruction error=0.28656295459729747, regularized loss=0.04105916347376638 regularized loss variation=0.007348268562704244.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.00034696406986296346, 0.01683129437953885]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=750, reconstruction error=0.28679901329718843, regularized loss=0.04112683701412043 regularized loss variation=0.007223567607007812.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0002627420826076366, 0.013158350081549685]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=800, reconstruction error=0.28642780296807, regularized loss=0.04102044315655777 regularized loss variation=0.008701879878733719.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0005447118439624145, 0.012428505553802023]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=850, reconstruction error=0.28648972323614685, regularized loss=0.04103818075996201 regularized loss variation=0.007251520740621105.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.00018992189106784157, 0.008296605896230323]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=900, reconstruction error=0.28640062165060637, regularized loss=0.04101265804092689 regularized loss variation=0.007326390485960887.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0003350501286559161, 0.008204618873378174]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=950, reconstruction error=0.2866152000984985, regularized loss=0.04107413646375117 regularized loss variation=0.007136617614327371.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0002553487604254914, 0.00696844039853362]
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.8554475564141224e-05, 2.080595490587331e-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.23452515071154423, regularized loss=0.027501023158136267 regularized loss variation=0.9533005344936518.
Feasibility gaps for A: [0.006934761900188146]
Feasibility gaps for the Bi-matrices: [0.47777479376893245, 0.17115970603542538]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=50, reconstruction error=0.10750045365430454, regularized loss=0.005778173767940639 regularized loss variation=0.0013405716899259517.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.00016623857103931355, 7.246566317706688e-05]
Feasibility gaps for C: [2.6506291951755746e-05]
Coupled matrix factorization iteration=100, reconstruction error=0.10419576717356265, regularized loss=0.005428378948443637 regularized loss variation=0.0009630728406895118.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.00016750078160195135, 3.616635894388375e-05]
Feasibility gaps for C: [7.26821720176826e-06]
Coupled matrix factorization iteration=150, reconstruction error=0.10267973576902208, regularized loss=0.005271564068798096 regularized loss variation=0.0003645893170536265.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [6.921858448254536e-05, 1.2105775335169636e-05]
Feasibility gaps for C: [1.5759648020172511e-06]
Coupled matrix factorization iteration=200, reconstruction error=0.10228923375381686, regularized loss=0.005231543670971494 regularized loss variation=4.591154652192792e-05.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.9530079003914102e-05, 3.491440592097337e-06]
Feasibility gaps for C: [4.266161848262022e-07]
Coupled matrix factorization iteration=250, reconstruction error=0.10223459686380783, regularized loss=0.005225956397952653 regularized loss variation=9.476502798857178e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [7.029351231086184e-06, 1.4496807014970646e-06]
Feasibility gaps for C: [1.4913567747234767e-07]
Coupled matrix factorization iteration=300, reconstruction error=0.10221835570407603, regularized loss=0.005224296121422506 regularized loss variation=4.762303505585412e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [3.3075002800552137e-06, 9.540021707179383e-07]
Feasibility gaps for C: [7.222596021095361e-08]
Coupled matrix factorization iteration=350, reconstruction error=0.10221112744662741, regularized loss=0.005223557286955356 regularized loss variation=7.102901331316581e-07.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [2.6716708962427e-06, 2.0909868264513295e-06]
Feasibility gaps for C: [3.492870238675301e-08]
Coupled matrix factorization iteration=400, reconstruction error=0.10220912796373718, regularized loss=0.005223352919553801 regularized loss variation=4.943398452114507e-07.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [8.440672894402484e-07, 2.509213055013502e-07]
Feasibility gaps for C: [1.8819619884105736e-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.554301975036419e-07, 1.5299990758210787e-07]
Feasibility gaps for C: [1.1117932795035532e-08]
Coupled matrix factorization iteration=500, reconstruction error=0.10220771274812757, regularized loss=0.0052232082726018795 regularized loss variation=1.594305610243083e-07.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [2.6714753081240184e-07, 9.7234784362899e-08]
Feasibility gaps for C: [7.240558357744665e-09]
Coupled matrix factorization iteration=550, reconstruction error=0.10220737567652786, regularized loss=0.00522317382134145 regularized loss variation=1.1046501417077894e-07.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.61935070531499e-07, 6.257492873085524e-08]
Feasibility gaps for C: [4.887272817845601e-09]
Coupled matrix factorization iteration=600, reconstruction error=0.10220713227080787, regularized loss=0.005223148943511208 regularized loss variation=8.299662168721872e-08.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.0046617909035999e-07, 4.103683375954143e-08]
Feasibility gaps for C: [3.3582157792555675e-09]
Coupled matrix factorization iteration=650, reconstruction error=0.10220694385991158, regularized loss=0.005223129686591559 regularized loss variation=6.604997199514979e-08.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [6.375153340120065e-08, 2.7991341266182157e-08]
Feasibility gaps for C: [2.321081839201791e-09]
Coupled matrix factorization iteration=700, reconstruction error=0.10220679073658931, regularized loss=0.005223114036336479 regularized loss variation=5.4718740445731665e-08.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [4.182358678221191e-08, 2.0641483045328327e-08]
Feasibility gaps for C: [1.5987344910846468e-09]
Coupled matrix factorization iteration=750, reconstruction error=0.10220666200327765, regularized loss=0.005223100878926119 regularized loss variation=4.661772100767867e-08.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [2.932120898693914e-08, 1.7034041699256395e-08]
Feasibility gaps for C: [1.0895583191623573e-09]
Coupled matrix factorization iteration=800, reconstruction error=0.10220655119384468, regularized loss=0.005223089553469997 regularized loss variation=4.0499093681892926e-08.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [2.306177814477151e-08, 1.5569054751104444e-08]
Feasibility gaps for C: [7.321781318522617e-10]
Coupled matrix factorization iteration=850, reconstruction error=0.10220645423101597, regularized loss=0.005223079643238381 regularized loss variation=3.566734138370049e-08.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [2.060625006473167e-08, 1.507346314235623e-08]
Feasibility gaps for C: [4.892420160889344e-10]
Coupled matrix factorization iteration=900, reconstruction error=0.1022063684021423, regularized loss=0.0052230708709772155 regularized loss variation=3.171457836942774e-08.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.998528026223582e-08, 1.4894580880469949e-08]
Feasibility gaps for C: [3.3951835166332577e-10]
Coupled matrix factorization iteration=950, reconstruction error=0.10220629181217035, regularized loss=0.00522306304299726 regularized loss variation=2.839062699402323e-08.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.9987943308000534e-08, 1.4749036103453516e-08]
Feasibility gaps for C: [2.689686256694736e-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.28228485865454467, regularized loss=0.039842370712808134 regularized loss variation=0.890970613788606.
Feasibility gaps for A: [0.009954018187787964]
Feasibility gaps for the Bi-matrices: [0.2239975451369923, 0.17952115245757755]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=50, reconstruction error=0.08754354867414268, regularized loss=0.0038319364572309946 regularized loss variation=0.002461984736989006.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0003200852385374236, 0.00029937909984316293]
Feasibility gaps for C: [5.4764882617851975e-06]
Coupled matrix factorization iteration=100, reconstruction error=0.08610152443153316, regularized loss=0.00370673625471695 regularized loss variation=0.000194677288040196.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [2.3242625553100167e-05, 7.56833562425614e-06]
Feasibility gaps for C: [6.01849186120492e-07]
Coupled matrix factorization iteration=150, reconstruction error=0.08566923630679649, regularized loss=0.003669609024694869 regularized loss variation=6.046423204695785e-05.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.3912219842697495e-05, 7.739185573320904e-06]
Feasibility gaps for C: [7.603533563821788e-07]
Coupled matrix factorization iteration=200, reconstruction error=0.08558989432012258, regularized loss=0.0036628150048648757 regularized loss variation=2.208073218952435e-05.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.703416411173623e-05, 1.372598713277024e-05]
Feasibility gaps for C: [4.5853278685136373e-07]
Coupled matrix factorization iteration=250, reconstruction error=0.08556226678342443, regularized loss=0.0036604507485589483 regularized loss variation=8.993080909365738e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [4.833122744559216e-06, 1.6707618331443756e-06]
Feasibility gaps for C: [2.529015985357008e-07]
Coupled matrix factorization iteration=300, reconstruction error=0.08554938728821843, regularized loss=0.0036593488326947944 regularized loss variation=4.391291929579049e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [3.4858066418549256e-06, 9.816368913181838e-07]
Feasibility gaps for C: [1.6875814714739563e-07]
Coupled matrix factorization iteration=350, reconstruction error=0.0855422090171864, regularized loss=0.003658734761770003 regularized loss variation=2.5749473006509268e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [2.400633896824215e-06, 7.05169999028699e-07]
Feasibility gaps for C: [1.1048268571381145e-07]
Coupled matrix factorization iteration=400, reconstruction error=0.08553776232221537, regularized loss=0.0036583543915459037 regularized loss variation=1.7472727775609412e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.6895476517012008e-06, 4.5937765586484885e-07]
Feasibility gaps for C: [7.422848413890557e-08]
Coupled matrix factorization iteration=450, reconstruction error=0.08553458711919705, regularized loss=0.003658082796825755 regularized loss variation=1.32279478708762e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.1921430892069157e-06, 3.1996225250492426e-07]
Feasibility gaps for C: [4.857526549956372e-08]
Coupled matrix factorization iteration=500, reconstruction error=0.08553197246357841, regularized loss=0.0036578591567551677 regularized loss variation=1.1491760415282908e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [8.714455208716444e-07, 2.3382418693078474e-07]
Feasibility gaps for C: [3.391130211716336e-08]
Coupled matrix factorization iteration=550, reconstruction error=0.08553016663316992, regularized loss=0.003657704702148907 regularized loss variation=4.670535440578359e-07.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [6.80162156147631e-07, 3.145673538921224e-07]
Feasibility gaps for C: [2.7088409924488746e-08]
Coupled matrix factorization iteration=600, reconstruction error=0.08552942886193414, regularized loss=0.0036576416007243263 regularized loss variation=2.5468352524109695e-07.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [5.390204209685728e-07, 2.249831947037369e-07]
Feasibility gaps for C: [2.439864848988518e-08]
Coupled matrix factorization iteration=650, reconstruction error=0.08552901056564477, regularized loss=0.003657605824169087 regularized loss variation=1.4986421712216737e-07.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [3.9539260934777493e-07, 1.6535261819560842e-07]
Feasibility gaps for C: [1.9334167115712334e-08]
Coupled matrix factorization iteration=700, reconstruction error=0.08552875744488725, regularized loss=0.003657584175033178 regularized loss variation=9.367074595657998e-08.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [2.868526260664209e-07, 1.2268586363674203e-07]
Feasibility gaps for C: [1.4624891037846352e-08]
Coupled matrix factorization iteration=750, reconstruction error=0.08552859399001696, regularized loss=0.0036575701949545826 regularized loss variation=6.288414721265368e-08.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [2.0708048318517191e-07, 9.129176800732398e-08]
Feasibility gaps for C: [1.0785290468486386e-08]
Coupled matrix factorization iteration=800, reconstruction error=0.0855284793631205, regularized loss=0.0036575603910838643 regularized loss variation=4.661349634627937e-08.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.4500872920810291e-07, 6.506575458741322e-08]
Feasibility gaps for C: [7.723668105408273e-09]
Coupled matrix factorization iteration=850, reconstruction error=0.0855283888283115, regularized loss=0.0036575526477834195 regularized loss variation=3.924806165671953e-08.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.0071497002585228e-07, 4.686010931143251e-08]
Feasibility gaps for C: [5.348101445077418e-09]
Coupled matrix factorization iteration=900, reconstruction error=0.0855283086022975, regularized loss=0.0036575457861849183 regularized loss variation=3.6361579681482846e-08.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [7.172104615021364e-08, 3.511909056431399e-08]
Feasibility gaps for C: [3.6561134290231885e-09]
Coupled matrix factorization iteration=950, reconstruction error=0.0855282318618967, regularized loss=0.003657539222711181 regularized loss variation=3.567234405293391e-08.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [5.365912276447771e-08, 2.828551015714203e-08]
Feasibility gaps for C: [2.491141621907468e-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.38539255102332226, 0.13910021179771864]
Feasibility gaps for C: [0.0004440734510898378]
Coupled matrix factorization iteration=50, reconstruction error=0.10314779647046571, regularized loss=0.005319733958356309 regularized loss variation=0.016393928497468516.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.005516295343159788, 0.016255368366386636]
Feasibility gaps for C: [6.719203953685178e-06]
Coupled matrix factorization iteration=100, reconstruction error=0.08843959760502258, regularized loss=0.003910781212269158 regularized loss variation=0.00014158123803620687.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [4.021911960800775e-05, 1.522915030896427e-05]
Feasibility gaps for C: [2.1058715485317656e-07]
Coupled matrix factorization iteration=150, reconstruction error=0.08824914414838528, regularized loss=0.0038939557214612416 regularized loss variation=7.174068957244777e-05.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.52987757504635e-05, 6.136441967562542e-06]
Feasibility gaps for C: [4.8563022189210435e-08]
Coupled matrix factorization iteration=200, reconstruction error=0.08816767468728776, regularized loss=0.003886769429881701 regularized loss variation=9.307394454178884e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [3.1206811263566623e-06, 1.2561851492380753e-06]
Feasibility gaps for C: [7.232388055469572e-09]
Coupled matrix factorization iteration=250, reconstruction error=0.08815442354116529, regularized loss=0.0038856011949375785 regularized loss variation=4.303970118569062e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.5976449208981578e-06, 7.452299889144234e-07]
Feasibility gaps for C: [1.5649069582113144e-09]
Coupled matrix factorization iteration=300, reconstruction error=0.08814257038057799, regularized loss=0.003884556356647572 regularized loss variation=2.1065866156290174e-05.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [4.5243710822063925e-06, 1.7681147128027372e-06]
Feasibility gaps for C: [2.6334901289131926e-08]
Coupled matrix factorization iteration=350, reconstruction error=0.0881091994541177, regularized loss=0.0038816155142227474 regularized loss variation=1.1350540013463679e-05.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [3.7347852606215012e-06, 1.136449944883997e-06]
Feasibility gaps for C: [2.150812626345948e-08]
Coupled matrix factorization iteration=400, reconstruction error=0.08808925929962963, regularized loss=0.0038798588019786926 regularized loss variation=7.570921231673541e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.4464836332328609e-06, 4.624272095630052e-07]
Feasibility gaps for C: [5.281461271794322e-09]
Coupled matrix factorization iteration=450, reconstruction error=0.0880733428295357, regularized loss=0.003878456858584464 regularized loss variation=7.193449383834907e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [8.123825549273877e-07, 3.5684527727347305e-07]
Feasibility gaps for C: [3.937306714718955e-09]
Coupled matrix factorization iteration=500, reconstruction error=0.08805672089217481, regularized loss=0.003876993047141188 regularized loss variation=7.936801629061167e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [9.483883331608416e-07, 3.7384838986486083e-07]
Feasibility gaps for C: [6.454721479575378e-09]
Coupled matrix factorization iteration=550, reconstruction error=0.08803801817061611, regularized loss=0.0038753463217048667 regularized loss variation=9.08957337673315e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [9.945415176589418e-07, 3.4860163996414706e-07]
Feasibility gaps for C: [7.309879330855187e-09]
Coupled matrix factorization iteration=600, reconstruction error=0.08801646782983286, regularized loss=0.0038734493046200013 regularized loss variation=1.0493273620150191e-05.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.0750839306569611e-06, 3.463510403619336e-07]
Feasibility gaps for C: [7.656780689862051e-09]
Coupled matrix factorization iteration=650, reconstruction error=0.08799181974037792, regularized loss=0.0038712801706115812 regularized loss variation=1.1887851758096596e-05.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.1447032996978019e-06, 3.535997102802888e-07]
Feasibility gaps for C: [7.758898329611206e-09]
Coupled matrix factorization iteration=700, reconstruction error=0.0879647920856407, regularized loss=0.0038689023233349986 regularized loss variation=9.049174006836977e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.4523115152116663e-06, 7.434176208177889e-07]
Feasibility gaps for C: [7.346668141504085e-09]
Coupled matrix factorization iteration=750, reconstruction error=0.08794459286038715, regularized loss=0.0038671257066896295 regularized loss variation=9.001320229025535e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.119832971751489e-06, 3.2336679101396235e-07]
Feasibility gaps for C: [5.005284882487557e-09]
Coupled matrix factorization iteration=800, reconstruction error=0.08793516958658398, regularized loss=0.003866297025110642 regularized loss variation=8.162130417268884e-07.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [4.661512935396946e-07, 1.312870387968263e-07]
Feasibility gaps for C: [4.057845108610151e-09]
Coupled matrix factorization iteration=850, reconstruction error=0.0879335445216034, regularized loss=0.003866154126066404 regularized loss variation=6.824339617062457e-07.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [2.114626530071277e-07, 7.474506921715168e-08]
Feasibility gaps for C: [1.9456154199742023e-09]
Coupled matrix factorization iteration=900, reconstruction error=0.08793215525298956, regularized loss=0.0038660319637179295 regularized loss variation=5.859876222100645e-07.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.4762792153495584e-07, 5.245213325395404e-08]
Feasibility gaps for C: [1.21714028962037e-09]
Coupled matrix factorization iteration=950, reconstruction error=0.08793096213595326, regularized loss=0.003865927051077223 regularized loss variation=5.032383571702018e-07.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.0308682397354564e-07, 3.847838153520366e-08]
Feasibility gaps for C: [8.220238005870621e-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.2345704105542823, regularized loss=0.027511638753802276 regularized loss variation=0.9155336821898681.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.4363161673356732, 0.22188606285488285]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=50, reconstruction error=0.08574361603418466, regularized loss=0.0036759838453088445 regularized loss variation=0.0021871892088128177.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0002050877612682179, 4.7171206070219956e-05]
Feasibility gaps for C: [8.779796213431633e-06]
Coupled matrix factorization iteration=100, reconstruction error=0.08406783885619686, regularized loss=0.003533700764975741 regularized loss variation=0.0002469036627616713.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [5.743008602321173e-05, 9.416472593413833e-06]
Feasibility gaps for C: [1.2022907086936324e-06]
Coupled matrix factorization iteration=150, reconstruction error=0.08379989981858157, regularized loss=0.0035112116048021543 regularized loss variation=5.9432968376638225e-05.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [2.1258716470825822e-05, 3.5457986204746655e-06]
Feasibility gaps for C: [4.108290291755574e-07]
Coupled matrix factorization iteration=200, reconstruction error=0.08372954054052467, regularized loss=0.0035053179795636823 regularized loss variation=1.7131386796160946e-05.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.043486395248947e-05, 1.9735137265040806e-06]
Feasibility gaps for C: [1.920441951461934e-07]
Coupled matrix factorization iteration=250, reconstruction error=0.08370913304440829, regularized loss=0.0035036094775232236 regularized loss variation=4.848770590692833e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [5.176852570470948e-06, 9.981002990647986e-07]
Feasibility gaps for C: [9.150305945262662e-08]
Coupled matrix factorization iteration=300, reconstruction error=0.08370336249794326, regularized loss=0.003503126446731047 regularized loss variation=1.430772467640239e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [2.737556451031611e-06, 5.1730190932143e-07]
Feasibility gaps for C: [4.634201020821863e-08]
Coupled matrix factorization iteration=350, reconstruction error=0.08370168485312357, regularized loss=0.003502986023625808 regularized loss variation=3.9202219701621487e-07.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.3173624986729903e-06, 2.523209038527842e-07]
Feasibility gaps for C: [2.2113205055769665e-08]
Coupled matrix factorization iteration=400, reconstruction error=0.08370121255617433, regularized loss=0.0035029464916869375 regularized loss variation=1.2292429564117875e-07.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [6.510098400930565e-07, 1.1798808038315162e-07]
Feasibility gaps for C: [1.0738626913728203e-08]
Coupled matrix factorization iteration=450, reconstruction error=0.08370105022731258, regularized loss=0.003502932904577552 regularized loss variation=4.832062386367695e-08.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [3.265945553910173e-07, 5.869173668965424e-08]
Feasibility gaps for C: [5.372873824600436e-09]
Coupled matrix factorization iteration=500, reconstruction error=0.08370097936218071, regularized loss=0.0035029269730941005 regularized loss variation=2.3904839821785356e-08.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.642276574562458e-07, 2.9790160463018638e-08]
Feasibility gaps for C: [2.7153055756231015e-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.346179210652011e-08, 1.5642396175007713e-08]
Feasibility gaps for C: [1.3920417013623224e-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.143986624174163e-08, 1.0143057417168136e-08]
Feasibility gap for C: [8.6587789253999e-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.975029013458259

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: (0 minutes 51.509 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