Main Content

Implement Hardware-Efficient Complex Partial-Systolic Q-less QR with Forgetting Factor

This example shows how to use the hardware-efficient Complex Partial-Systolic Q-less QR Decomposition with Forgetting Factor block.

Q-less QR Decomposition with Forgetting Factor

The Complex Partial-Systolic Q-less QR Decomposition with Forgetting Factor block implements the following recursion to compute the upper-triangular factor R of continuously streaming n-by-1 row vectors A(k,:) using forgetting factor $\alpha$. It's as if matrix A is infinitely tall. The forgetting factor in the range $0<\alpha<1$ keeps it from integrating without bound.

$$&#xA;\begin{array}{ccl}&#xA;R_0 &#38; =&#38; \mbox{zeros}(n,n)\\[2ex]&#xA;[\sim\;,\;R_1] &#38;=&#38; \mbox{qr}\left(\alpha\left[\begin{array}{c}R_0\\&#xA;A(1,:)\end{array}\right],\; 0\right)\\[4ex]&#xA;[\sim\;,\;R_2] &#38;=&#38; \mbox{qr}\left(\alpha\left[\begin{array}{c}R_1\\&#xA;A(2,:)\end{array}\right],\; 0\right)\\[4ex]&#xA;\vdots\\[4ex]&#xA;[\sim\;,\;R_k] &#38;=&#38; \mbox{qr}\left(\alpha\left[\begin{array}{c}R_{k-1}\\&#xA;A(k,:)\end{array}\right],\; 0\right)\\[4ex]&#xA;\vdots\\[4ex]&#xA;\end{array}&#xA;$$

Define System Parameters

n is the length of the row vectors A(k,:) and the number of rows and columns in R.

n = 5;

m is the effective numbers of rows of A to integrate over.

m = 100;

Use the fixed.forgettingFactor function to compute the forgetting factor as a function of the number of rows that you are integrating over.

forgettingFactor = fixed.forgettingFactor(m)
forgettingFactor =

    0.9950

precisionBits defines the number of bits of precision required for the QR Decomposition. Set this value according to system requirements.

precisionBits = 24;

In this example, complex-valued matrix A is constructed such that the magnitude of the real and imaginary parts of its elements is less than or equal to one, so the maximum possible absolute value of any element is $|1+1i|=\sqrt{2}$. Your own system requirements will define what those values are. If you don't know what they are, and A is a fixed-point input to the system, then you can use the upperbound function to determine the upper bounds of the fixed-point types of A.

max_abs_A is an upper bound on the maximum magnitude element of A.

max_abs_A = sqrt(2);

Select Fixed-Point Types

Use the fixed.qlessqrFixedpointTypes function to compute fixed-point types.

T = fixed.qlessqrFixedpointTypes(m,max_abs_A,precisionBits)
T = 

  struct with fields:

    A: [0x0 embedded.fi]

T.A is the fixed-point type computed for transforming A to R in-place so that it does not overflow.

T.A
ans = 

[]

          DataTypeMode: Fixed-point: binary point scaling
            Signedness: Signed
            WordLength: 30
        FractionLength: 24

Define Simulation Parameters

Create random matrix A to contain a specified number of inputs.

numInputs is the number of input rows A(k,:) for this example.

numInputs = 500;
rng('default')
A = fixed.example.complexUniformRandomArray(-1,1,numInputs,n);

Cast the inputs to the types determined by fixed.qlessqrFixedpointTypes.

A = cast(A,'like',T.A);

Cast the forgetting factor to a fixed-point type with the same word length as A and best-precision scaling.

forgettingFactor = fi(forgettingFactor,1,T.A.WordLength);

Select a stop time for the simulation that is long enough to process all the inputs from A.

stopTime = 4*numInputs*T.A.WordLength;

Open the Model

model = 'ComplexPartialSystolicQlessQRForgettingFactorModel';
open_system(model);

Set Variables in the Model Workspace

Use the helper function setModelWorkspace to add the variables defined above to the model workspace.

fixed.example.setModelWorkspace(model,'A',A,'n',n,...
    'forgettingFactor',forgettingFactor,...
    'stopTime',stopTime);

Simulate the Model

out = sim(model);

Verify the Accuracy of the Output

Define matrix $A_k$ as follows

$$A_k = \left[\begin{array}{cccc}\alpha^k \\&#38; \alpha^{k-1} \\&#xA;&#38; &#38; \ddots\\&#38;&#xA;&#38; &#38; \alpha\end{array}\right]A(1:k,\; :).$$

Then using the formula for the computation of the $k$ th output $R_k$, and the fact that $[Q,R]=\mbox{qr}(A,0) \Rightarrow A'A = R'Q'QR = R'R$, you can show that

$$A_k'A_k = R_k'R_k.$$

So to verify the output, the difference between $A_k'A_k$ and $R_k'R_k$ should be small.

Choose the last output of the simulation.

R = double(out.R(:,:,end))
R =

  Columns 1 through 4

   7.8315 + 0.0000i   0.2537 + 0.1037i  -0.1812 + 0.1835i   0.7054 + 0.2411i
   0.0000 + 0.0000i   7.7917 + 0.0000i  -0.2939 - 0.6734i  -0.2813 - 0.2824i
   0.0000 + 0.0000i   0.0000 + 0.0000i   8.2110 + 0.0000i  -0.2445 - 0.0531i
   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   8.0730 + 0.0000i
   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i

  Column 5

  -0.2351 - 0.1305i
  -0.4660 - 0.6805i
  -0.2956 + 0.5699i
   0.2053 - 0.0133i
   8.0363 + 0.0000i

Verify that R is upper triangular.

isequal(R,triu(R))
ans =

  logical

   1

Verify that the diagonal is greater than or equal to zero.

diag(R)
ans =

    7.8315
    7.7917
    8.2110
    8.0730
    8.0363

Synchronize the last output R with the input by finding the number of inputs that produced it.

A = double(A);
alpha = double(forgettingFactor);
relative_errors = nan(1,n);
for k = 1:numInputs
    A_k = alpha.^(k:-1:1)' .* A(1:k,:);
    relative_errors(k) = norm(A_k'*A_k - R'*R)/norm(A_k'*A_k);
end

k is the number of inputs A(k,:) that produced the last R.

k = find(relative_errors==min(relative_errors),1,'last')
k =

   497

Verify that

$$A_k'A_k = R_k'R_k$$

with a small relative error.

A_k =  alpha.^(k:-1:1)' .* A(1:k,:);
relative_error = norm(A_k'*A_k - R'*R)/norm(A_k'*A_k)
relative_error =

   6.3862e-06

Suppress mlint warnings in this file.

%#ok<*NOPTS>

See Also