Rank-One Nonnegative Matrix Factorization

Author

CVXPY Developers and Balasubramanian Narasimhan

Introduction

The DGP atom library has several functions of positive matrices, including the trace, (matrix) product, sum, Perron-Frobenius eigenvalue, and $(I - X)^{- 1}$ (eye-minus-inverse). In this example, we use some of these atoms to approximate a partially known elementwise positive matrix as the outer product of two positive vectors.

We would like to approximate $A$ as the outer product of two positive vectors $x$ and $y$ , with $x$ normalized so that the product of its entries equals $1$ . Our criterion is the average relative deviation between the entries of $A$ and $x y^{T}$ , that is,

$\frac{1}{m n} \sum_{i = 1}^{m} \sum_{j = 1}^{n} R (A_{i j}, x_{i} y_{j}),$

where $R$ is the relative deviation of two positive numbers, defined as

$R (a, b) = max {a / b, b / a} - 1.$

The corresponding optimization problem is

$\begin{array}{ll} minimize & \frac{1}{m n} \sum_{i = 1}^{m} \sum_{j = 1}^{n} R (X_{i j}, x_{i} y_{j}) \\ subject to & x_{1} x_{2} \dots x_{m} = 1 \\ X_{i j} = A_{i j}, for (i, j) \in Ω, \end{array}$

with variables $X \in R_{+ +}^{m \times n}$ , $x \in R_{+ +}^{m}$ , and $y \in R_{+ +}^{n}$ . We can cast this problem as an equivalent generalized geometric program by discarding the $- 1$ from the relative deviations.

Problem Data

We consider the partially known matrix

$A = [\begin{matrix} 1.0 & ? & 1.9 \\ ? & 0.8 & ? \\ 3.2 & 5.9 & ? \end{matrix}],$

where the question marks denote the missing entries.

Problem Formulation and Solution

m <- 3
n <- 3

X <- Variable(c(m, n), pos = TRUE)
x <- Variable(m, pos = TRUE)
y <- Variable(n, pos = TRUE)

# Construct the outer product matrix: each row i is x[i] * y^T
outer_product <- vstack(t(x[1] * y), t(x[2] * y), t(x[3] * y))

# Relative deviations (without the -1): max(X/outer_product, outer_product/X)
relative_deviations <- max_elemwise(
  multiply(X, power(outer_product, -1)),
  multiply(power(X, -1), outer_product)
)

Warning: `multiply()` is deprecated. Use `x * y` instead.
This warning is displayed once per session.

objective <- sum_entries(relative_deviations)

constraints <- list(
  X[1, 1] == 1.0,
  X[1, 3] == 1.9,
  X[2, 2] == 0.8,
  X[3, 1] == 3.2,
  X[3, 2] == 5.9,
  x[1] * x[2] * x[3] == 1.0
)

problem <- Problem(Minimize(objective), constraints)
result <- psolve(problem, gp = TRUE)
check_solver_status(problem)

cat(sprintf("Optimal value (average relative deviation): %g\n",
            1.0 / (m * n) * (result - m * n)))
cat("Outer product approximation:\n")
print(value(outer_product))
cat("x:", value(x), "\n")
cat("y:", value(y), "\n")

Optimal value (average relative deviation): 0
Outer product approximation:
          [,1]    [,2]      [,3]
[1,] 1.0000000 1.84375 1.9000000
[2,] 0.4338983 0.80000 0.8244068
[3,] 3.2000000 5.90000 6.0800000
x: 0.8963701 0.3889335 2.868384 
y: 1.115611 2.056907 2.11966

The optimal average relative deviation is essentially zero, meaning the outer product $x y^{T}$ matches $A$ exactly at the known entries. The solver fills in the unknown entries of $X$ to be consistent with the rank-one outer product structure.

Session Info

R version 4.5.2 (2025-10-31)
Platform: aarch64-apple-darwin20
Running under: macOS Tahoe 26.3

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: America/Los_Angeles
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] CVXR_1.8.1

loaded via a namespace (and not attached):
 [1] slam_0.1-55       cli_3.6.5         knitr_1.51        ECOSolveR_0.6.1  
 [5] rlang_1.1.7       xfun_0.56         clarabel_0.11.2   otel_0.2.0       
 [9] gurobi_13.0-1     Rglpk_0.6-5.1     highs_1.12.0-3    cccp_0.3-3       
[13] scs_3.2.7         S7_0.2.1          jsonlite_2.0.0    Rcplex_0.3-8     
[17] backports_1.5.0   rprojroot_2.1.1   htmltools_0.5.9   Rmosek_11.1.1    
[21] gmp_0.7-5.1       piqp_0.6.2        rmarkdown_2.30    grid_4.5.2       
[25] evaluate_1.0.5    fastmap_1.2.0     yaml_2.3.12       compiler_4.5.2   
[29] codetools_0.2-20  htmlwidgets_1.6.4 Rcpp_1.1.1        here_1.0.2       
[33] osqp_1.0.0        lattice_0.22-9    digest_0.6.39     checkmate_2.3.4  
[37] Matrix_1.7-4      tools_4.5.2

References

Agrawal, A., Diamond, S., Boyd, S. (2019). Disciplined Geometric Programming. Optimization Letters, 13(5), 961–976.