Network topology from Spectral templates code always results in NaN for lambda

[thildebrant]

I am trying to use matlab and YALMIP to solve a graph learning problem of recovering eigenvalues from the eigenvectors of the covariance of sampled graph signal data. This is to implement the algorithm in https://graphex.mit.edu/sites/default/files/documents/2017_Poster_Segarra.pdf. My objective is to estimate the Adjacency matrix.
I have tried different versions of the constraints and other troubleshooting, but lambda is always NaN or very small, and estimated_A is also very small, so interpreted as 0.
This is my first time working through this kind of optimization problem and I appreciate any assistance. Thank you.

My code is

% Generate a random network with 100 nodes
N = 100;
G = gsp_sensor(N);
A = G.A;
D = diag(sum(A, 2));
L = D - A;   % Laplacian matrix

% Generate observed graph signals
X = sin(rand(N, N));

% Diffuse graph signals using the heat equation
alpha = .5;
X_diff = zeros(N, N);
for t = 1:N
    X_diff(:, t) = (eye(N) - alpha * L) * X(:, t);
end

% Compute sample covariance matrix
C = (X_diff * X_diff') / N;

% Extract spectral templates by computing the eigenvectors and eigenvalues of the sample covariance matrix
[V, D] = eig(C);

% Formulate optimization problem with YALMIP
A_est = sdpvar(N, N, 'full');  % Unknown graph shift operator (Laplacian matrix)
lambda = sdpvar(N, 1);       % Unknown eigenvalues

options = sdpsettings('verbose', 2, 'solver', 'sedumi');
% Define the objective function (minimize L1 norm of A)
objective = norm(A_est(:), 1);

% Set up the constraints
constraints = [A_est == V * diag(lambda) * V', A_est == A_est', diag(A_est) == zeros(N, 1), A_est >= 0, sum(A_est, 2) ==1];

% Define objective function
objective = norm(A_est(:), 1);

% Solve SDP problem using YALMIP
sol = solvesdp(constraints, objective, options);

disp(value(A_est));
disp(value(lambda));

[johanlofberg]

First, it is bit unclear what you actually look for. You've explicitly said A_est should have a non-symmetric parameterization, but then added the constraint that it is symmetric, and then A_est>=0. I hope you realize that means you are looking for a non-negative symmetric matrix (i.e. not positive semidefinite). Regardless, it is a bit strange to force a use of twice as many variables as required. If you want a symmetric non-negative matrix, you can just use A_est=sdpvar(N) and then impose the non-negativity with, e.g., triu(A_est)>=0

The only way for lambda to contain NaN is that there are zero columns/rows in V meaning some lamda aren't part of of the problem (i.e. any value is just as valid). A_est being 0 simply cannot happen if the problem has been solved, as you have the constraint that the rows sum to 1.

[johanlofberg]

and there is no official command in YALMIP called solvesdp https://yalmip.github.io/tutorial/basics/

[johanlofberg]

and when I solve it for a generic psd matrix A, the problem is generically infeasible

The fact that you tell YALMIP to use sedumi, despite the fact that you've set up a linear program, indicates that you are not modelling what you actually want

[thildebrant]

Thank you for the feedback. Let me try to address your points:
A_est should be symmetric, as the graph I am looking for should be undirected and unweighted, ideally with no self-loops (diagonal = 0). I am trying to find weights (lambda) that minimize the norm of A_est, given that I know the eigenvectors from the sampled data.
Concerning the solver: you are correct, I used sedumi because of "Linear programming can be solved by quadratic programming which can be solved by second-order cone programming which can be solved by semidefinite programming". Based on the feedback would you have any suggestions? (The paper mentioned that it's possible to determine the graph laplacian in a similar manner just with different constraint, but I was looking at the adjacency matrix first as I thought it was simpler)

[johanlofberg]

Hence A_est = sdpvar(N) with triu(A_est)>=0

There must be something structurally wrong with the model as it is infeasible (unless you data looks very particular compared to the random data I generated). Also, you haven't told us that the solver actually has returned infeasibility status on your model, but that is my guess

Feels overdetermined, With a symmetric parameterization, A_est will have N(N+1)/2 parameters, and lambda has N, but you have N(N+1)/2+N+1 equality constraints, i.e. more equality constraints than parameters