% from ChatGPT, this works pretty well

clear
clc

A = double(imread('peppers.bmp'));
[m,n] = size(A);

k = 100;  % <-- MUCH BETTER for 256×256

% Rank-k SVD
[U,S,V] = svds(A, k);

% Leverage scores
row_lev = sum(U.^2,2);
col_lev = sum(U.^2,2);

row_prob = row_lev / sum(row_lev);
col_prob = col_lev / sum(col_lev);

% row_idx = randsample(m, k, true, row_prob);
% col_idx = randsample(n, k, true, col_prob);

row_idx = randsample(m, k, true, row_prob);
col_idx = randsample(n, k, true, col_prob);

C = A(:, col_idx);
R = A(row_idx, :);
W = A(row_idx, col_idx);

U_cur = pinv(W);
A_cur = C * U_cur * R;

imshow(uint8(A_cur));

myPSNR = psnr(A, A_cur, 255);
mySSIM = ssim(A, A_cur);