SVD to a matrix subset (rows)

Hello everyone

I have a quite large matrix A (MxN with M>N) and a for cycle of many (tens of thousands) iterations where, for each iteration, I have to exclude some rows of the matrix A and perform a single value decomposition of it.

Something like

for n = 1:NLoops
    [U,S,V] = svd( A(rows_selection(n,:),:) );
    
    % Do stuff with U, S and V...
end

where rows_selection is a logical matrix that select the rows to be analyzed for each loop.

Since this code happens to be really slow, is it possible to calculate [U,S,V] = svd(A) before the for cycle and then arrange the U, S and V matrices based on the particular rows to be selected for each for loop?

Answer 1

[John D'Errico]

No, I don't think you can easily find the SVD of a subset of the rows of a matrix, based on the SVD of the full matrix. At least, not easily. (@Christine Tobler) may have an idea.

But whenever I see a question like this, I wonder if the real problem is, you are trying to solve the wrong problem. That is, one reason for computing the SVD repeatedly may be to find the numerical rank of each subset of the rows. Or maybe you want to find the subset of rows that are maximally independent from each other. And if you are doing something like one of these things, there are likely better ways of doing so than repeat calls to SVD.

It is often the case that someone asks a question, looking to make the scheme they have run faster, but the real issue is, you need to improve the approach you are taking. Find a better way of solving the problem than what you are doing.

Edit:

For example, if your goal is to find the subset of k rows of your matrix, such that the condition number of the resulting matrix is as small as possible, when you extract that set of rows and then compute the condition number on that subset? Now, rather than computing cond on EVERY subset of rows, I might try using GA. Have GA select the subset, with minimal condition number. Will this really gain? Well, yes. Consider this example:

A = rand(1e6,10);

Now, my goal is to find the subset of 10 rows from all 1e6 rows, such that the condition number is as small as possible. It might be a moderately interesting problem (though I have no idea if it is your real problem.) Were we to apply brute force, that would necessitate MATLAB to compute

nchoosek(1e6,10)

Warning: Result may not be exact. Coefficient has a maximum relative error of 3.5527e-15, corresponding to absolute error 9.789885939637849e+38.

ans = 2.7556e+53

calls to COND (which calls SVD). And clearly that would be an insane thing to do. But we could use GA instead, to identify a subset of those rows, with a minimal condition number. Yes, GA might not find the perfectly optimal solution, but it would be massively better than anything you could do using brute force, and the result would take only some thousands of calls to COND. Even if it was many thousands of calls, you win by use of an optimization tool here.

Agin, I have no clue as to why you have decided to perform this programming task. Maybe someone told you to do it this way, or maybe you decided it was the way to solve some problem. But that does not make necessarily a good way to approach the real underlying problem. We don't know, because all we see is the final question that you think to be important.

[John D'Errico]

Just for kicks, I decided it might be interesting to write code to solve the problem I posed. But since I'm a lazy guy, why write code when a computer can do it for me? So I posed the question to the MATLAB chat AI tool. Here is what it wrote. (Yes, it used globals. One serious demerit for that, when there was no need to do so. And there were a couple of bugs in the code it wrote. And it never checks to see if there were any duplicates in the row selection. That would generate an infinite condition number, so a poor choice of rows. More demerits. But overall, the code will almost run, after I fixed the multiple bugs in the code it wrote.)

% Define the matrix A
% Global variable for the matrix A
global A;
A = rand(100, 5); % Example: 100 rows and 5 columns

% Number of rows to select
k = 5;

% Objective function to minimize the condition number
function condNum = conditionNumberObjective(selectedRows)
    global A;
    subMatrix = A(selectedRows,:);
    condNum = cond(subMatrix); % Calculate the condition number
end

% Set up the genetic algorithm
nvars = k; % Number of variables (rows to select)
IntCon = 1:nvars; % All variables are integers
lb = ones(1, k); % Lower bounds
ub = size(A, 1) * ones(1, k); % Upper bounds

% Ensure unique selections
options = optimoptions('ga', 'PopulationSize', 100, 'MaxGenerations', 200, 'Display', 'iter', 'UseParallel', false, 'CreationFcn', @gacreationuniform);


A = A; % Assign the matrix to the global variable

% Run genetic algorithm
[selectedRows, fval] = ga(@conditionNumberObjective, nvars, [], [], [], [], lb, ub, [], IntCon);

ga stopped because the average change in the penalty function value is less than options.FunctionTolerance and 
the constraint violation is less than options.ConstraintTolerance.

% Display the results
disp('Selected Rows:');

Selected Rows:

disp(selectedRows);

    83    34    66    25    74

disp('Minimum Condition Number:');

Minimum Condition Number:

disp(fval);

    4.7195

Now, had I just chosen some arbitrary subset of 5 rows, what might I see for a condition number?

ind = randi(100,[1,5])

ind = 1×5

75 24 15 5 53

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

cond(A(ind,:))

ans = 22.8638

Anyway, I still have no clue as to why you want to do what you are doing. But I will conjecture you may be doing it for the wrong reason.

[Davide Cataldo]

Thank you for your time John.

What I am trying to do is implementing an algorithm for radar interferometry application (phase unwrapping algorithm from "Monserrat, O.. (2012). Deformation measurement and monitoring with GB-SAR." if you are interested).

In short,

1) I perform a least square estimation (the paper proposes SVD-LS estimation) of N variables from M>N observations.

2) Calculate the residuals between estimations and observations.

3) Perform SVD-LS estimation again without considering those observations for which the residuals exceed a threshold (thus, removing some rows from the matrix A).

4) Calculate the residuals again.

5) Evaluate which observations (those that were excluded at step 3) are to be definitively excluded, can be corrected or reintegrated.

6) Repeat untill all (or enough) residuals remain within the threshold.

7) Repeat for all the pixels of interest in the image under test.

TLDR: I have to perform a least square estimation using the entire matrix A first and then, iteratively, remove rows from the matrix (and perform LS estimation again) untill I find the optimal set of observations (subset of rows).

Another option was to use the lsqminnorm Matlab function, which is slightly faster than the LS estimation using the SVD.

[Steven Lord]

What are the values of M and N for the problem you're considering? Knowing how big a problem you're trying to solve may help us more easily assess if any approaches we are about to suggest are feasible; possible but going to take a very, very long time; or not possible in the lifetime of this universe.

Answer 2

[Chuguang Pan]

I find there is a svdappend function for calculating SVD of matrix incrementally, such as from SVD of A to SVD of [A,D] or [A;D], without explicitly forming full matrix [A,D] or [A;D]. However, the for loop of your requirement is not appending rows, maybe the logical selection of different rows can be transformed into appending form so that the svdappend function is useful.

[Torsten]

The existence of "svdappend" seems to indicate that there exist at least techniques for "updating" or "downdating" an SVD when adding/removing columns. Thus a completely new SVD doesn't seem necessary in this case.