Implement Hardware-Efficient Complex Burst QR Decomposition

This example demonstrates how to compute the QR decomposition of complex matrices using hardware-efficient MATLAB® code in Simulink®. This model shares computational resources across steps of the QR Decomposition algorithm. It thus uses fewer on chip resources than a fully pipelined approach, while sacrificing some total throughput.

Solving Matrix Equations with the QR Decomposition

When solving matrix equations, it is seldom, if ever, necessary to compute the inverse of a matrix [1][3]. The Complex Burst QR Decomposition block provides a hardware efficient way of solving the equation

where is a complex m x n matrix, is a complex n x p matrix, and is a complex m x p matrix. This equation can be converted to the form

through a series of orthogonal transformations. Here, is a complex m x n upper triangular matrix. is a complex m x m orthogonal matrix such that . The Complex Burst QR Decomposition retains only the non-zero rows of , along with the corresponding rows of .

The QR Decomposition is well suited to fixed-point architectures because it can be entirely performed with Givens rotations. These have an efficient fixed-point implementation in terms of the CORDIC Algorithm . For more details on the algorithm for fixed-point QR decomposition, see Perform QR Factorization Using CORDIC. The Simulink model used in this example is:

fxpdemo_complex_burst_QR_model

This model works with fixed-point, double, single, and scaled double data types.

Input Parameters for Complex Burst QR Decomposition

The Complex Burst QR Decomposition block takes four input parameters. The mask for this block is shown below:

You must assign input data to variables in the MATLAB workspace. To customize this algorithm, change the data according to the procedures used in the following sections.

Define Matrix Dimensions

The model's architecture allocates the minimal memory necessary to store the data needed to perform the QR decomposition. Therefore, the size of the input matrices must be known at compile time. Here, m is the number of rows in the complex matrices A and B, n is the number of columns in A, and p is the number of columns in B.

Additionally, since the block can handle decomposition of an indefinite number of matrices in series, a number of samples is specified. The model terminates when all of the samples have been used.

m = 4;
n = 4;
p = 4;
num_samples = 100;

Define Input Datatypes

Before generating input data, it is important to specify the data type of the matrix data as shown below. This example uses signed fixed-point types with 16 bit word lengths and fractional types. For a complex 4 x 4 matrix, the output type needs an additional 3 bits to accommodate the data growth in the integer part of the input data; see Use Hardware-Efficient Algorithm to Solve Systems of Complex-Valued Linear Equations for further details.

word_length = 16;
qr_growth_bits = 3;
fraction_length = word_length - 1;
nt = numerictype(1,word_length + qr_growth_bits,fraction_length)

Define the Number of CORDIC Iterations

The CORDIC approximation for Givens rotations is an iterative algorithm that gives an additional bit of accuracy per iteration. For signed fixed-point datatypes, the greatest possible accuracy is achieved when the number of iterations performed is one fewer than the wordlength.

NumberOfCORDICIterations = nt.WordLength - 1

Initialize Random Data

The data handler in this example takes complex matrices A and B as inputs. In this example, A and B are defined as random matrices elements drawn from a uniform distribution from -1 to 1. Note that A and B are each defined as three-dimensional arrays, where each sample matrix is stored in the first two dimensions.

A = fi(complex(2*rand(m,n,num_samples) - 1, 2*rand(m,n,num_samples) - 1),nt);
B = fi(complex(2*rand(m,p,num_samples) - 1, 2*rand(m,p,num_samples) - 1),nt);

Transferring Data to Complex Burst QR Decomposition

The ready port triggers the Data Handler subsystem. When ready is high, the block asserts validIn and sends the next row of A and B to aIn and bIn. The protocol for transferring data allows data to be sent whenever ready is high, ensuring all data is processed. If data is sent when ready is not high, it will not be processed.

Handling Outputs

Complex Burst QR Decomposition outputs data one row at a time. Whenever the block outputs a result row, it asserts validOut. Note that rOut outputs the rows of , while cOut outputs the rows of . The rows are output in reverse order, as this is the natural order to consume them for back substitution. In this example, they are saved as (m * num_samples) x n matrices, with the rows ordered as they were received.

Simulate

sim fxpdemo_complex_burst_QR_model

Reconstruct the Solutions from the Output Data

Because the data are output in reverse order, you must reconstruct the data to interpret the results. The following code puts the data in the correct order.

C = cell(1,num_samples);
R = cell(1,num_samples);
for ii = 1:num_samples
    C{ii} = flipud(C_Out((ii-1)*m + 1:ii*m,:));
    R{ii} = flipud(R_Out((ii-1)*m + 1:ii*m,:));
end

Evaluate the Accuracy of the Result

To evaluate the accuracy of the Burst QR Decomposition, examine the magnitude of the difference between the solution to the matrix equation using the 'Complex QR Burst Decomposition' block, and that obtained using MATLAB®'s built-in backsolve for doubles. Plot the absolute error and condition number for each sample. The data show that the accuracy of the solution tracks the condition number, as expected [2].

xAct = cell(1,num_samples);
xExp = cell(1,num_samples);
xErr = zeros(1,num_samples);
condNumber = zeros(1,num_samples);
for ii = 1:num_samples
    xAct{ii} = double(R{ii})\double(C{ii});
    xExp{ii} = double(A(:,:,ii))\double(B(:,:,ii));
    xErr(ii) = norm(xAct{ii} - xExp{ii});
    condNumber(ii) = cond(double(A(:,:,ii)));
end
figure(1)
clf
h1 = subplot(2,1,1);
hold on;
h1.YScale = 'log';
plot(xErr)
grid on
title('Error of ''Complex Burst QR Decomposition''');
ylabel('norm(R\\C - A\\B)');
h2 = subplot(2,1,2);
hold on;
h2.YScale = 'log';
plot(condNumber);
grid on
title('Condition Number of Samples');
ylabel('cond(A)');
xlabel('Test Point');
linkaxes([h1,h2],'x');

Close the model

close_system fxpdemo_complex_burst_QR_model

References

[1] George E. Forsythe, M.A. Malcom and Cleve B. Moler. Computer Methods for Mathematical Computations. Englewood Cliffs, N.J.: Prentice-Hall, 1977.

[2] Cleve B. Moler. Cleve's Corner: What is the Condition Number of a Matrix?, The MathWorks, Inc. 2017.

[3] Cleve B. Moler. Numerical Computing with MATLAB. SIAM, 2004. isbn: 978-0-898716-60-3.

%#ok<*NASGU,*NOPTS>