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.3314035749797264, regularized loss=0.05491416475467156 regularized loss variation=0.9067503005970162.
Feasibility gaps for A: [0.007604248629946112]
Feasibility gaps for the Bi-matrices: [0.49058444441575316, 0.21620083723066494]
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.005599748831554872, 0.07667375386621973]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=100, reconstruction error=0.29785648881423793, regularized loss=0.04435924396437312 regularized loss variation=0.0011881407930269376.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.005009444420332494, 0.07180217320018313]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=150, reconstruction error=0.294138000996816, regularized loss=0.04325858181520146 regularized loss variation=0.0036106456523246305.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.003470803899617987, 0.07184365518506443]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=200, reconstruction error=0.2906407963526251, regularized loss=0.04223603625224404 regularized loss variation=0.0064873535427517354.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0012741643853616715, 0.04740122551753363]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=250, reconstruction error=0.2894397507992609, regularized loss=0.04188768467136912 regularized loss variation=0.006209991163836618.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0006535342383363337, 0.029785117592718284]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=300, reconstruction error=0.2889316012131703, regularized loss=0.041740735089803235 regularized loss variation=0.004945541949489507.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0006393959965885707, 0.023528987537592916]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=350, reconstruction error=0.28829963601881436, regularized loss=0.04155834006429042 regularized loss variation=0.005982919305416159.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0004594696378956256, 0.016428340281644527]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=400, reconstruction error=0.2883317467091478, regularized loss=0.04156759808017408 regularized loss variation=0.0033198151656166942.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.00032441720970438734, 0.011851148272540458]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=450, reconstruction error=0.2883962258944582, regularized loss=0.041586191555083675 regularized loss variation=0.005117950303896945.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0003424531287140414, 0.010356397804950337]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=500, reconstruction error=0.28801236795810914, regularized loss=0.041475562048418624 regularized loss variation=0.0044134307782602725.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0001994316444438696, 0.007295906159875317]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=550, reconstruction error=0.28745181764455224, regularized loss=0.04131427373357846 regularized loss variation=0.009270091648538247.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.00023002428324835433, 0.006620289322964216]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=600, reconstruction error=0.28817402162912936, regularized loss=0.04152213337095296 regularized loss variation=0.0032028417183609975.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.00015090801583757948, 0.005146935692336602]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=650, reconstruction error=0.28820775591846937, regularized loss=0.04153185528578001 regularized loss variation=0.005204518314978592.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.00016502926783164767, 0.004844781998492579]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=700, reconstruction error=0.2876713406601516, regularized loss=0.041377400118604496 regularized loss variation=0.005534934380316159.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0001236100808969609, 0.007050310433727524]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=750, reconstruction error=0.2883352098128162, regularized loss=0.04156859660890037 regularized loss variation=0.0006546385849540619.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0001320590418935247, 0.004003453585765779]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=800, reconstruction error=0.28811365302017755, regularized loss=0.041504738528315634 regularized loss variation=0.003340502622031597.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.00011390540491884459, 0.0032805137953514653]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=850, reconstruction error=0.288049114498451, regularized loss=0.04148614618167087 regularized loss variation=0.0060685182771680135.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.00010613331409364367, 0.003102418995939569]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=900, reconstruction error=0.28714333912288625, regularized loss=0.04122564860132043 regularized loss variation=0.007963031513012444.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0005503475581184724, 0.0027588118695728034]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=950, reconstruction error=0.2882115699720609, regularized loss=0.04153295453288007 regularized loss variation=0.0012041085479763867.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.00017316401850567929, 0.002637321649711234]
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.00999730042828625]
Feasibility gaps for the Bi-matrices: [0.20663569486247202, 0.21316798877616955]
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.0008513823031412011, 0.0008772327964805434]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=100, reconstruction error=0.30225051736801, regularized loss=0.04567768762461486 regularized loss variation=6.670798920721441e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [6.668298424645227e-05, 5.083663270520364e-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.4967275268850378e-05, 1.4205180473683634e-05]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=200, reconstruction error=0.3021764997407128, regularized loss=0.0456553184977745 regularized loss variation=1.1125916035544285e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.512459639877436e-05, 3.1009754457098596e-05]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=250, reconstruction error=0.3021186234043592, regularized loss=0.04563783130387251 regularized loss variation=7.136904699411424e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [5.202263077518944e-05, 4.56335953343618e-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.38228852648131817, 0.23073647629374147]
Feasibility gaps for C: [0.0009314483755838897]
Coupled matrix factorization iteration=50, reconstruction error=0.3123683026075554, regularized loss=0.048786978236962646 regularized loss variation=0.006813974457742268.
Feasibility gaps for A: [1.062493357055203e-05]
Feasibility gaps for the Bi-matrices: [0.004992585559639324, 0.055370055664303516]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=100, reconstruction error=0.31193360650529645, regularized loss=0.04865128743370056 regularized loss variation=0.005804161385191129.
Feasibility gaps for A: [8.270516437516441e-06]
Feasibility gaps for the Bi-matrices: [0.0034403594488344115, 0.047167333037834214]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=150, reconstruction error=0.31182068055500745, regularized loss=0.048616068410894 regularized loss variation=0.00560518384506813.
Feasibility gaps for A: [6.239341621301804e-06]
Feasibility gaps for the Bi-matrices: [0.002548771440675913, 0.04327512747482581]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=200, reconstruction error=0.3116668434349995, regularized loss=0.048568110648368246 regularized loss variation=0.005913499828937766.
Feasibility gaps for A: [4.7529316400836964e-06]
Feasibility gaps for the Bi-matrices: [0.0019970316023714347, 0.04107598244417706]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=250, reconstruction error=0.3116204925213335, regularized loss=0.04855366567961923 regularized loss variation=0.005018590608084654.
Feasibility gaps for A: [3.7733997928012866e-06]
Feasibility gaps for the Bi-matrices: [0.0017128224561445903, 0.040541435113762236]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=300, reconstruction error=0.31138518454562175, regularized loss=0.048480366577255456 regularized loss variation=0.003815024716722892.
Feasibility gaps for A: [3.3675121451255395e-06]
Feasibility gaps for the Bi-matrices: [0.00139849160972112, 0.03928419291081779]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=350, reconstruction error=0.31132090414112956, regularized loss=0.04846035267762519 regularized loss variation=0.005458365600897585.
Feasibility gaps for A: [2.1920479712771665e-06]
Feasibility gaps for the Bi-matrices: [0.0012215913812588597, 0.03905084679401711]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=400, reconstruction error=0.3110550421867215, regularized loss=0.04837761963489154 regularized loss variation=0.002270918878000745.
Feasibility gaps for A: [2.1463021214521014e-06]
Feasibility gaps for the Bi-matrices: [0.0010863750845886524, 0.03833256770073943]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=450, reconstruction error=0.3111549612164016, regularized loss=0.0484087049447902 regularized loss variation=0.0031276563852613132.
Feasibility gaps for A: [1.78811124794108e-06]
Feasibility gaps for the Bi-matrices: [0.001041032512737332, 0.038262180695189586]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=500, reconstruction error=0.3109198357089825, regularized loss=0.04833557211865034 regularized loss variation=0.002581446382635471.
Feasibility gaps for A: [1.5045328206120823e-06]
Feasibility gaps for the Bi-matrices: [0.000976412654976418, 0.03794814576805481]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=550, reconstruction error=0.311069373226882, regularized loss=0.0483820774798826 regularized loss variation=0.003321014993842006.
Feasibility gaps for A: [1.3022301154719827e-06]
Feasibility gaps for the Bi-matrices: [0.0009567921989914468, 0.0379309467314847]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=600, reconstruction error=0.31085595799316484, regularized loss=0.04831571330992413 regularized loss variation=0.002705006294939972.
Feasibility gaps for A: [1.1320479560719347e-06]
Feasibility gaps for the Bi-matrices: [0.0009217250949397318, 0.03776161765288152]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=650, reconstruction error=0.3110232586340207, regularized loss=0.04836773370566246 regularized loss variation=0.0034104324022427752.
Feasibility gaps for A: [1.0047060161538166e-06]
Feasibility gaps for the Bi-matrices: [0.0009125304803454128, 0.037759971581629914]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=700, reconstruction error=0.31081834077137405, regularized loss=0.048304020479935 regularized loss variation=0.002768274543813954.
Feasibility gaps for A: [8.936164345928955e-07]
Feasibility gaps for the Bi-matrices: [0.0008916682259105455, 0.03765781410847938]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=750, reconstruction error=0.3109944071458426, regularized loss=0.04835876063799705 regularized loss variation=0.0034600559241850515.
Feasibility gaps for A: [8.07533041004766e-07]
Feasibility gaps for the Bi-matrices: [0.0008871280615395349, 0.037661001439576236]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=800, reconstruction error=0.3107935465899643, regularized loss=0.04829631430098415 regularized loss variation=0.0028056609032922645.
Feasibility gaps for A: [7.303411821203237e-07]
Feasibility gaps for the Bi-matrices: [0.0008737596804719195, 0.037594066443221975]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=850, reconstruction error=0.3109747012239086, regularized loss=0.0483526324006496 regularized loss variation=0.0034908177567477897.
Feasibility gaps for A: [6.689459559521199e-07]
Feasibility gaps for the Bi-matrices: [0.0008714571877195043, 0.03759858154582901]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=900, reconstruction error=0.31077602180761266, regularized loss=0.04829086786528287 regularized loss variation=0.002829910193987209.
Feasibility gaps for A: [6.126679939502112e-07]
Feasibility gaps for the Bi-matrices: [0.0008623633343147749, 0.037551936098245094]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=950, reconstruction error=0.3109604239955354, regularized loss=0.048348192645741576 regularized loss variation=0.0035114397779978505.
Feasibility gaps for A: [5.670366508391938e-07]
Feasibility gaps for the Bi-matrices: [0.0008612050870965695, 0.037556561633287154]
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.008775850160562646]
Feasibility gaps for the Bi-matrices: [0.4583428412993907, 0.2812577162876004]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=50, reconstruction error=0.2939515799183767, regularized loss=0.04320376566825491 regularized loss variation=0.004044211601716679.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0016129311723681696, 0.05593480922354094]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=100, reconstruction error=0.29340842031307723, regularized loss=0.043044250555307695 regularized loss variation=0.004083989387752174.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.001864057664228692, 0.05100739577781305]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=150, reconstruction error=0.2931404510511532, regularized loss=0.04296566202123678 regularized loss variation=0.003998429724350731.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0012148515386552984, 0.046771554489357665]
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.001138353345621646, 0.04603530584659063]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=250, reconstruction error=0.29302641978196015, regularized loss=0.042932241345116764 regularized loss variation=0.0038705031359493112.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0010761947515037795, 0.045071978106245886]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=300, reconstruction error=0.2929899080343091, regularized loss=0.042921543104976444 regularized loss variation=0.0038114040386139185.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.001013365524566322, 0.04387127525024219]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=350, reconstruction error=0.29296417215027337, regularized loss=0.04291400308184751 regularized loss variation=0.003762958445477458.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0009525989176797745, 0.04245084717816407]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=400, reconstruction error=0.29294909977861205, regularized loss=0.0429095875305496 regularized loss variation=0.0037247879957865476.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0008934570527588717, 0.040839461516833465]
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.0014827875191874894, 0.045925674961990126]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=500, reconstruction error=0.2885811626607134, regularized loss=0.04163954372130456 regularized loss variation=0.008388168588552.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0025613957285298542, 0.06203337425244079]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=550, reconstruction error=0.28743614027640096, regularized loss=0.041309767368497424 regularized loss variation=0.006051730404067943.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0014052168420058237, 0.04158232217837158]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=600, reconstruction error=0.28675716124665807, regularized loss=0.04111483476312093 regularized loss variation=0.0074778029575014344.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0008192613861116502, 0.031227676788273793]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=650, reconstruction error=0.28691263481554086, regularized loss=0.04115943000839795 regularized loss variation=0.007220804238638106.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0005198054696550911, 0.02314656180833286]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=700, reconstruction error=0.2865629545972972, regularized loss=0.04105916347376631 regularized loss variation=0.007348268562705973.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.00034696406986298503, 0.016831294379538748]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=750, reconstruction error=0.28679901329718793, regularized loss=0.04112683701412029 regularized loss variation=0.0072235676070076665.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0002627420826076712, 0.01315835008154979]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=800, reconstruction error=0.286427802968069, regularized loss=0.041020443156557486 regularized loss variation=0.008701879878725173.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0005447118439624238, 0.0124285055538021]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=850, reconstruction error=0.28648972323614735, regularized loss=0.041038180759962156 regularized loss variation=0.007251520740619193.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.00018992189106782327, 0.008296605896230323]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=900, reconstruction error=0.28640062165060615, regularized loss=0.041012658040926825 regularized loss variation=0.007326390485961241.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.00033505012865591804, 0.008204618873378035]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=950, reconstruction error=0.2866152000984983, regularized loss=0.041074136463751104 regularized loss variation=0.00713661761433338.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.00025534876042550286, 0.006968440398533555]
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.8554475564005258e-05, 2.0805954905903676e-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.340422185442039

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.23452515071154362, regularized loss=0.027501023158136125 regularized loss variation=0.9533005344936523.
Feasibility gaps for A: [0.006934761900188113]
Feasibility gaps for the Bi-matrices: [0.47777479376893267, 0.1711597060354255]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=50, reconstruction error=0.1075004536543052, regularized loss=0.00577817376794071 regularized loss variation=0.0013405716899388101.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0001662385710393658, 7.246566317703268e-05]
Feasibility gaps for C: [2.6506291951779396e-05]
Coupled matrix factorization iteration=100, reconstruction error=0.10419576717355984, regularized loss=0.005428378948443346 regularized loss variation=0.0009630728407563835.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.0001675007816019134, 3.616635894387809e-05]
Feasibility gaps for C: [7.2682172017337696e-06]
Coupled matrix factorization iteration=150, reconstruction error=0.10267973576902277, regularized loss=0.005271564068798167 regularized loss variation=0.00036458931705362155.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [6.921858448260267e-05, 1.2105775335254953e-05]
Feasibility gaps for C: [1.575964801998636e-06]
Coupled matrix factorization iteration=200, reconstruction error=0.1022892337538133, regularized loss=0.0052315436709711284 regularized loss variation=4.591154659172425e-05.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.9530079003895627e-05, 3.4914405919081893e-06]
Feasibility gaps for C: [4.26616184821548e-07]
Coupled matrix factorization iteration=250, reconstruction error=0.10223459686380783, regularized loss=0.005225956397952653 regularized loss variation=9.476502784749837e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [7.029351231253105e-06, 1.4496807015423243e-06]
Feasibility gaps for C: [1.49135677486414e-07]
Coupled matrix factorization iteration=300, reconstruction error=0.10221835570407317, regularized loss=0.005224296121422214 regularized loss variation=4.7623035756473975e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [3.3075002801459525e-06, 9.540021706153389e-07]
Feasibility gaps for C: [7.222596021413881e-08]
Coupled matrix factorization iteration=350, reconstruction error=0.1022111274466267, regularized loss=0.0052235572869552835 regularized loss variation=7.102901470796883e-07.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [2.671670896175618e-06, 2.090986826406109e-06]
Feasibility gaps for C: [3.4928702400113455e-08]
Coupled matrix factorization iteration=400, reconstruction error=0.10220912796373575, regularized loss=0.005223352919553655 regularized loss variation=4.943398594921453e-07.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [8.440672893235596e-07, 2.509213055148703e-07]
Feasibility gaps for C: [1.8819619870036996e-08]
Coupled matrix factorization iteration=450, reconstruction error=0.10220822705576495, regularized loss=0.005223260838941401 regularized loss variation=2.572463903227659e-07.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [4.5543019759947314e-07, 1.5299990769358177e-07]
Feasibility gaps for C: [1.1117932814648264e-08]
Coupled matrix factorization iteration=500, reconstruction error=0.10220771274812829, regularized loss=0.005223208272601953 regularized loss variation=1.594305047302584e-07.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [2.6714753088150877e-07, 9.723478440971512e-08]
Feasibility gaps for C: [7.240558381321315e-09]
Coupled matrix factorization iteration=550, reconstruction error=0.10220737567652714, regularized loss=0.005223173821341376 regularized loss variation=1.1046502828590169e-07.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.6193507038622152e-07, 6.257492880686705e-08]
Feasibility gaps for C: [4.887272822766922e-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.0046617907011548e-07, 4.103683377286898e-08]
Feasibility gaps for C: [3.3582157559800277e-09]
Coupled matrix factorization iteration=650, reconstruction error=0.10220694385991014, regularized loss=0.005223129686591411 regularized loss variation=6.605005618841688e-08.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [6.375153338725305e-08, 2.7991341226444377e-08]
Feasibility gaps for C: [2.321081827744945e-09]
Coupled matrix factorization iteration=700, reconstruction error=0.10220679073659003, regularized loss=0.005223114036336553 regularized loss variation=5.471871221516225e-08.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [4.1823586644546034e-08, 2.0641482978432382e-08]
Feasibility gaps for C: [1.5987344977011618e-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.9321208849162032e-08, 1.7034041764228345e-08]
Feasibility gaps for C: [1.0895583148126208e-09]
Coupled matrix factorization iteration=800, reconstruction error=0.10220655119384899, regularized loss=0.005223089553470436 regularized loss variation=4.04989674740463e-08.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [2.30617781920416e-08, 1.556905473701787e-08]
Feasibility gaps for C: [7.321781146544963e-10]
Coupled matrix factorization iteration=850, reconstruction error=0.10220645423101454, regularized loss=0.0052230796432382355 regularized loss variation=3.56673973470229e-08.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [2.060625008374367e-08, 1.507346308975525e-08]
Feasibility gaps for C: [4.892420419573385e-10]
Coupled matrix factorization iteration=900, reconstruction error=0.10220636840214087, regularized loss=0.00522307087097707 regularized loss variation=3.171459248483025e-08.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.998528041361966e-08, 1.48945808363587e-08]
Feasibility gaps for C: [3.395183447005287e-10]
Coupled matrix factorization iteration=950, reconstruction error=0.10220629181216893, regularized loss=0.005223063042997115 regularized loss variation=2.839068295752411e-08.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.9987943290333445e-08, 1.4749036086283782e-08]
Feasibility gaps for C: [2.68968624890026e-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.2822848586545441, regularized loss=0.039842370712807974 regularized loss variation=0.8909706137886064.
Feasibility gaps for A: [0.009954018187787999]
Feasibility gaps for the Bi-matrices: [0.22399754513699158, 0.17952115245757755]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=50, reconstruction error=0.08754354867414019, regularized loss=0.003831936457230776 regularized loss variation=0.002461984737102778.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.00032008523853818623, 0.00029937909984414273]
Feasibility gaps for C: [5.476488261806104e-06]
Coupled matrix factorization iteration=100, reconstruction error=0.08610152443153485, regularized loss=0.0037067362547170963 regularized loss variation=0.00019467728796124538.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [2.324262555292627e-05, 7.568335624290742e-06]
Feasibility gaps for C: [6.018491861253156e-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.3912219842721752e-05, 7.739185573313899e-06]
Feasibility gaps for C: [7.603533563826365e-07]
Coupled matrix factorization iteration=200, reconstruction error=0.08558989432012687, regularized loss=0.003662815004865243 regularized loss variation=2.208073204945997e-05.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.7034164111704996e-05, 1.3725987132724086e-05]
Feasibility gaps for C: [4.585327868783503e-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.833122744548068e-06, 1.6707618331757488e-06]
Feasibility gaps for C: [2.52901598560243e-07]
Coupled matrix factorization iteration=300, reconstruction error=0.08554938728821843, regularized loss=0.0036593488326947944 regularized loss variation=4.39129188928494e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [3.4858066417752897e-06, 9.816368913625648e-07]
Feasibility gaps for C: [1.6875814716456655e-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.4006338968399806e-06, 7.051699989723414e-07]
Feasibility gaps for C: [1.1048268572997162e-07]
Coupled matrix factorization iteration=400, reconstruction error=0.08553776232221451, regularized loss=0.00365835439154583 regularized loss variation=1.7472728972913525e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.6895476515843272e-06, 4.5937765594337526e-07]
Feasibility gaps for C: [7.422848413073443e-08]
Coupled matrix factorization iteration=450, reconstruction error=0.08553458711919705, regularized loss=0.003658082796825755 regularized loss variation=1.3227948672300082e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.1921430890721879e-06, 3.199622525739903e-07]
Feasibility gaps for C: [4.857526547225756e-08]
Coupled matrix factorization iteration=500, reconstruction error=0.0855319724635767, regularized loss=0.0036578591567550216 regularized loss variation=1.1491761009274569e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [8.714455207766742e-07, 2.3382418691961563e-07]
Feasibility gaps for C: [3.391130211357286e-08]
Coupled matrix factorization iteration=550, reconstruction error=0.08553016663316823, regularized loss=0.003657704702148762 regularized loss variation=4.6705354346502297e-07.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [6.801621562868614e-07, 3.145673539625475e-07]
Feasibility gaps for C: [2.7088409932777148e-08]
Coupled matrix factorization iteration=600, reconstruction error=0.08552942886193328, regularized loss=0.0036576416007242526 regularized loss variation=2.546835253596705e-07.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [5.39020420969504e-07, 2.2498319479212942e-07]
Feasibility gaps for C: [2.439864847244993e-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.9539260916324864e-07, 1.6535261813604506e-07]
Feasibility gaps for C: [1.933416711787205e-08]
Coupled matrix factorization iteration=700, reconstruction error=0.0855287574448864, regularized loss=0.003657584175033106 regularized loss variation=9.367080524172291e-08.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [2.868526260757085e-07, 1.2268586362638173e-07]
Feasibility gaps for C: [1.4624891049256546e-08]
Coupled matrix factorization iteration=750, reconstruction error=0.08552859399001782, regularized loss=0.003657570194954656 regularized loss variation=6.28841671325372e-08.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [2.0708048302669854e-07, 9.129176801500345e-08]
Feasibility gaps for C: [1.0785290503849944e-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.4500872907846322e-07, 6.506575453881002e-08]
Feasibility gaps for C: [7.72366809982605e-09]
Coupled matrix factorization iteration=850, reconstruction error=0.08552838882831235, regularized loss=0.003657552647783492 regularized loss variation=3.9248061893861393e-08.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.0071497002010787e-07, 4.686010929190826e-08]
Feasibility gaps for C: [5.3481014439647374e-09]
Coupled matrix factorization iteration=900, reconstruction error=0.08552830860229409, regularized loss=0.0036575457861846264 regularized loss variation=3.6361659480132347e-08.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [7.172104615655237e-08, 3.5119090567535064e-08]
Feasibility gaps for C: [3.6561134560533813e-09]
Coupled matrix factorization iteration=950, reconstruction error=0.08552823186189841, regularized loss=0.0036575392227113273 regularized loss variation=3.5672284292768866e-08.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [5.36591227350829e-08, 2.82855101599485e-08]
Feasibility gaps for C: [2.491141661231208e-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.22477174868817332, regularized loss=0.02526116950416967 regularized loss variation=0.9300354886525044.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.38539255102332054, 0.13910021179771861]
Feasibility gaps for C: [0.0004440734510898233]
Coupled matrix factorization iteration=50, reconstruction error=0.103147796470465, regularized loss=0.005319733958356236 regularized loss variation=0.01639392849746874.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.005516295343159798, 0.01625536836638708]
Feasibility gaps for C: [6.719203953677072e-06]
Coupled matrix factorization iteration=100, reconstruction error=0.08843959760502175, regularized loss=0.003910781212269085 regularized loss variation=0.00014158123801780642.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [4.021911960789612e-05, 1.5229150308966595e-05]
Feasibility gaps for C: [2.1058715480328947e-07]
Coupled matrix factorization iteration=150, reconstruction error=0.08824914414838611, regularized loss=0.0038939557214613153 regularized loss variation=7.174068955351574e-05.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.529877575039149e-05, 6.13644196754705e-06]
Feasibility gaps for C: [4.856302218711016e-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.1206811262422514e-06, 1.2561851493809306e-06]
Feasibility gaps for C: [7.23238807460607e-09]
Coupled matrix factorization iteration=250, reconstruction error=0.08815442354116694, regularized loss=0.003885601194937724 regularized loss variation=4.303970062205177e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.5976449210493663e-06, 7.452299889345629e-07]
Feasibility gaps for C: [1.5649069431496225e-09]
Coupled matrix factorization iteration=300, reconstruction error=0.0881425703805813, regularized loss=0.0038845563566478644 regularized loss variation=2.1065866099911593e-05.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [4.524371082311128e-06, 1.768114712726675e-06]
Feasibility gaps for C: [2.633490130000159e-08]
Coupled matrix factorization iteration=350, reconstruction error=0.08810919945412102, regularized loss=0.0038816155142230397 regularized loss variation=1.1350539976035161e-05.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [3.734785260687703e-06, 1.1364499447945194e-06]
Feasibility gaps for C: [2.150812627475581e-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.4464836333524827e-06, 4.624272095979399e-07]
Feasibility gaps for C: [5.281461263780513e-09]
Coupled matrix factorization iteration=450, reconstruction error=0.08807334282953819, regularized loss=0.0038784568585846827 regularized loss variation=7.19344938416995e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [8.123825550903418e-07, 3.5684527730752896e-07]
Feasibility gaps for C: [3.937306724987937e-09]
Coupled matrix factorization iteration=500, reconstruction error=0.08805672089217315, regularized loss=0.003876993047141041 regularized loss variation=7.936801629173327e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [9.483883331448895e-07, 3.738483897888382e-07]
Feasibility gaps for C: [6.454721489065473e-09]
Coupled matrix factorization iteration=550, reconstruction error=0.08803801817061528, regularized loss=0.0038753463217047934 regularized loss variation=9.089573414557418e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [9.94541517700635e-07, 3.4860163990442335e-07]
Feasibility gaps for C: [7.309879313240081e-09]
Coupled matrix factorization iteration=600, reconstruction error=0.08801646782983037, regularized loss=0.0038734493046197823 regularized loss variation=1.0493273657769377e-05.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.0750839305158007e-06, 3.463510402988009e-07]
Feasibility gaps for C: [7.656780680943254e-09]
Coupled matrix factorization iteration=650, reconstruction error=0.08799181974037958, regularized loss=0.0038712801706117265 regularized loss variation=1.1887851720568607e-05.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.144703299617414e-06, 3.535997102388774e-07]
Feasibility gaps for C: [7.7588983363363e-09]
Coupled matrix factorization iteration=700, reconstruction error=0.08796479208564155, regularized loss=0.003868902323335073 regularized loss variation=9.049173950006157e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.4523115153535518e-06, 7.434176207877064e-07]
Feasibility gaps for C: [7.34666813942658e-09]
Coupled matrix factorization iteration=750, reconstruction error=0.08794459286038797, regularized loss=0.0038671257066897015 regularized loss variation=9.001320286106414e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.1198329719210347e-06, 3.233667909317646e-07]
Feasibility gaps for C: [5.005284882428269e-09]
Coupled matrix factorization iteration=800, reconstruction error=0.0879351695865873, regularized loss=0.003866297025110934 regularized loss variation=8.162128529457512e-07.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [4.6615129365402845e-07, 1.3128703880482718e-07]
Feasibility gaps for C: [4.057845098334453e-09]
Coupled matrix factorization iteration=850, reconstruction error=0.08793354452160257, regularized loss=0.0038661541260663307 regularized loss variation=6.824339996209224e-07.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [2.1146265310006826e-07, 7.474506920492806e-08]
Feasibility gaps for C: [1.9456154037133006e-09]
Coupled matrix factorization iteration=900, reconstruction error=0.08793215525298956, regularized loss=0.0038660319637179295 regularized loss variation=5.859876032521298e-07.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.4762792163238229e-07, 5.245213318797799e-08]
Feasibility gaps for C: [1.2171403239435207e-09]
Coupled matrix factorization iteration=950, reconstruction error=0.08793096213595243, regularized loss=0.0038659270510771502 regularized loss variation=5.032384332283804e-07.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.0308682409976265e-07, 3.847838160165407e-08]
Feasibility gaps for C: [8.220238232479966e-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.23457041055428293, regularized loss=0.027511638753802425 regularized loss variation=0.9155336821898676.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.43631616733567374, 0.22188606285488283]
Feasibility gaps for C: [0.0]
Coupled matrix factorization iteration=50, reconstruction error=0.08574361603418551, regularized loss=0.0036759838453089173 regularized loss variation=0.0021871892088127743.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [0.00020508776126829377, 4.717120607029569e-05]
Feasibility gaps for C: [8.779796213409317e-06]
Coupled matrix factorization iteration=100, reconstruction error=0.08406783885619686, regularized loss=0.003533700764975741 regularized loss variation=0.00024690366280288726.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [5.743008602340581e-05, 9.416472593421727e-06]
Feasibility gaps for C: [1.202290708715281e-06]
Coupled matrix factorization iteration=150, reconstruction error=0.0837998998185807, regularized loss=0.003511211604802081 regularized loss variation=5.943296843962369e-05.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [2.1258716470613814e-05, 3.5457986204513547e-06]
Feasibility gaps for C: [4.1082902914041356e-07]
Coupled matrix factorization iteration=200, reconstruction error=0.08372954054052381, regularized loss=0.0035053179795636103 regularized loss variation=1.7131386795913867e-05.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.0434863952617311e-05, 1.9735137265034673e-06]
Feasibility gaps for C: [1.9204419513817687e-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.17685257056471e-06, 9.981002990475454e-07]
Feasibility gaps for C: [9.150305946242782e-08]
Coupled matrix factorization iteration=300, reconstruction error=0.08370336249794501, regularized loss=0.0035031264467311933 regularized loss variation=1.4307724885620213e-06.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [2.737556451040671e-06, 5.173019093308462e-07]
Feasibility gaps for C: [4.634201020066314e-08]
Coupled matrix factorization iteration=350, reconstruction error=0.0837016848531227, regularized loss=0.0035029860236257346 regularized loss variation=3.920221970162231e-07.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.3173624989426956e-06, 2.523209037744228e-07]
Feasibility gaps for C: [2.2113205063944694e-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.510098400991339e-07, 1.1798808039937233e-07]
Feasibility gaps for C: [1.0738626920764286e-08]
Coupled matrix factorization iteration=450, reconstruction error=0.08370105022731258, regularized loss=0.003502932904577552 regularized loss variation=4.832058238898102e-08.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [3.265945555128549e-07, 5.869173675501091e-08]
Feasibility gaps for C: [5.372873850552464e-09]
Coupled matrix factorization iteration=500, reconstruction error=0.08370097936218245, regularized loss=0.0035029269730942458 regularized loss variation=2.3904777423924174e-08.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [1.6422765769652605e-07, 2.979016056160081e-08]
Feasibility gaps for C: [2.7153055668933174e-09]
Coupled matrix factorization iteration=550, reconstruction error=0.08370094131155674, regularized loss=0.003502923788220333 regularized loss variation=1.3929728159435487e-08.
Feasibility gaps for A: [0.0]
Feasibility gaps for the Bi-matrices: [8.346179221336527e-08, 1.5642396161922354e-08]
Feasibility gaps for C: [1.3920417120980392e-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.143986623475755e-08, 1.0143057502946161e-08]
Feasibility gap for C: [8.658778620972456e-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.9750290134582587

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 51.873 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