How to solve a linear system for only one component in MATLAB

Question

I need to solve the linear system

A x = b

which can be done efficiently by

x = A \ b

But now A is very large and I actually only need one component, say x(1). Is there a way to solve this more efficiently than to compute all components of x?

A is not sparse. Here, efficiency is actually an issue because this is done for many b.

Also, storing the inverse of K and multiplying only its first row to b is not possible because K is badly conditioned. Using the \ operator employs the LDL solver in this case, and accuracy is lost when the inverse is explicitly used.

score 2 · Answer 1 · 2018-05-31T00:47:30.873

I don't think you'd technically get a speed-up over the very optimized Matlab routine however if you understand how it is solved then you can just solve for one part of x. E.g the following. in traditional solver you use backsub for QR solve for instance. In LU solve you use both back sub and front sub. I could get LU. Unfortunately, it actually starts at the end due to how it solves it. The same is true for LDL which would employ both. That doesn't preclude that fact there may be more efficient ways of solving whatever you have.

 function [Q,R] = qrcgs(A)
%Classical Gram Schmidt for an m x n matrix 

[m,n] = size(A); 
% Generates the Q, R matrices
Q = zeros(m,n); 
R = zeros(n,n);
    for k = 1:n
        % Assign the vector for normalization
        w = A(:,k);
        for j=1:k-1
            % Gets R entries
            R(j,k) = Q(:,j)'*w;
        end

        for j = 1:k-1
           % Subtracts off orthogonal projections 
           w = w-R(j,k)*Q(:,j);
        end
        % Normalize 
        R(k,k) = norm(w);
        Q(:,k) = w./R(k,k);
    end

end

function x = backsub(R,b)
% Backsub for upper triangular matrix.
[m,n] = size(R);
p = min(m,n);
x  = zeros(n,1);

    for i=p:-1:1
        % Look from bottom, assign to vector
        r = b(i);
        for j=(i+1):p
            % Subtract off the difference
            r = r-R(i,j)*x(j);
        end
        x(i) = r/R(i,i);
    end

end

That's a nice and basic answer. However I will not try it out because as you say it will most likely still be slower than the built-in methods. — rehctawrats, Jun 01 '18 at 10:44

OriolAbril · Accepted Answer · 2018-06-01T11:06:50.617

The method mldivide, generally represented as \ accepts solving many systems with the same A at once.

x = A\[b1 b2 b3 b4] # where bi are vectors with n rows

Solves the system for each b, and will return an nx4 matrix, where each column is the solution of each b. Calling mldivide like this should improve efficiency becaus the descomposition is only done once.

As in many decompositions like LU od LDL' (and in the one you are interested in particular) the matrix multiplying x is upper diagonal, the first value to be solved is x(n). However, having to do the LDL' decomposition, a simple backwards substitution algorithm won't be the bottleneck of the code. Therefore, the decomposition can be saved in order to avoid repeating the calculation for every bi. Thus, the code would look similar to this:

[LA,DA] = ldl(A);
DA = sparse(DA);
% LA = sparse(LA); %LA can also be converted to sparse matrix
% loop over bi
  xi = LA'\(DA\(LA\bi));
% end loop

As you can see in the documentation of mldivide (Algorithms section), it performs some checks on the input matrixes, and having defined LA as full and DA as sparse, it should directly go for a triangular solver and a tridiagonal solver. If LA was converted to sparse, it would use a triangular solver too, and I don't know if the conversion to sparse would represent any improvement.

Storing the factorization makes a lot of sense. I've actually done it before in another context but couldn't think of it now. Thanks for the reminder. — rehctawrats, Jun 01 '18 at 10:41
Happy to help, the key is that `\` uses many different solvers depending on the input matrix, as explained in the Algorithms section, so applying mldivide to the alredy decomposed matrixes should be much faster than on the original system. — OriolAbril, Jun 01 '18 at 11:04

How to solve a linear system for only one component in MATLAB

2 Answers2