fixed.jacobiSVD

Fixed-point Jacobi singular value decomposition

Since R2023a

Syntax

S = fixed.jacobiSVD(A)

[U,s,V] = fixed.jacobiSVD(A)

[___] = fixed.jacobiSVD(___,numberOfSweeps)

Description

S = fixed.jacobiSVD(A) returns a vector containing the singular values of matrix A in descending order.

[U,s,V] = fixed.jacobiSVD(A) performs a singular value decomposition of matrix A such that A = U*diag(s)*V'. s is a vector of nonnegative elements in decreasing order. U and V are unitary matrices.

[___] = fixed.jacobiSVD(___,numberOfSweeps) performs numberOfSweeps Jacobi iterations. If numberOfSweeps is not supplied, the default is 10.

Examples

Compute Singular Values of Fixed-Point Matrix

Compute the singular values of a full rank scaled-double matrix.

A = fi([1 0 1; -1 -2 0; 0 1 -1]);

Define fixed-point types that will never overflow. First, use the fixed.singularValueUpperBound function to determine the upper bound on the singular values. Then, define the integer length based on the value of the upper bound, with one additional bit for the sign, another additional bit for intermediate CORDIC growth, and one more bit for intermediate growth to compute the Jacobi rotations. Compute the fraction length based on the integer length and the desired word length.

svdUpperBound = fixed.singularValueUpperBound(3,3,max(abs(A(:))));
additionalBitGrowth = 3;
integerLength = ceil(log2(svdUpperBound)) + additionalBitGrowth;
wordLength = 16;
fractionLength = wordLength - integerLength;

Cast the matrix A to the signed fixed-point type.

T.A = fi([],1,wordLength,fractionLength,'DataType','Fixed');
A = cast(A,'like',T.A)

A = 
     1     0     1
    -1    -2     0
     0     1    -1

          DataTypeMode: Fixed-point: binary point scaling
            Signedness: Signed
            WordLength: 16
        FractionLength: 10

Compute the singular values of the fixed-point matrix A.

s = fixed.jacobiSVD(A)

s = 
    2.4482
    1.6865
    0.2295

          DataTypeMode: Fixed-point: binary point scaling
            Signedness: Signed
            WordLength: 16
        FractionLength: 10

The singular values are returned in a column vector in decreasing order with the same data type as the input matrix A.

Find Fixed-Point Singular Value Decomposition

Find the singular value decomposition of the fixed-point matrix A.

m = 4;
n = 2;
rng('default');
A = 10*randn(m,n)

A = 4×2

    5.3767    3.1877
   18.3389  -13.0769
  -22.5885   -4.3359
    8.6217    3.4262

Define fixed-point types that will never overflow. First, use the fixed.singularValueUpperBound function to determine the upper bound on the singular values.Then, define the integer length based on the value of the upper bound, with one additional bit for the sign, another additional bit for intermediate CORDIC growth, and one more bit for intermediate growth to compute the Jacobi rotations. Compute the fraction length based on the integer length and the desired word length.

svdUpperBound = fixed.singularValueUpperBound(m,n,max(abs(A(:))));
additionalBitGrowth = 3;
integerLength = ceil(log2(svdUpperBound)) + additionalBitGrowth;
wordLength = 16;
fractionLength = wordLength - integerLength;

Cast the matrix A to the signed fixed-point type.

T.A = fi([],1,wordLength,fractionLength,'DataType','Fixed');
A = cast(A,'like',T.A)

A = 
    5.3750    3.1875
   18.3359  -13.0781
  -22.5859   -4.3359
    8.6250    3.4297

          DataTypeMode: Fixed-point: binary point scaling
            Signedness: Signed
            WordLength: 16
        FractionLength: 7

Find the singular value decomposition of the fixed-point matrix A.

[U,s,V] = fixed.jacobiSVD(A)

U = 
   -0.1584    0.2739
   -0.6402   -0.7533
    0.7042   -0.5052
   -0.2618    0.3195

          DataTypeMode: Fixed-point: binary point scaling
            Signedness: Signed
            WordLength: 16
        FractionLength: 14

s = 
   31.0312
   14.1484

          DataTypeMode: Fixed-point: binary point scaling
            Signedness: Signed
            WordLength: 16
        FractionLength: 7

V = 
   -0.9918    0.1276
    0.1276    0.9918

          DataTypeMode: Fixed-point: binary point scaling
            Signedness: Signed
            WordLength: 16
        FractionLength: 14

Confirm the relation A = U*diag(s)*V'.

U*diag(s)*V'

ans = 
    5.3691    3.2170
   18.3441  -13.1050
  -22.5841   -4.3018
    8.6355    3.4472

          DataTypeMode: Fixed-point: binary point scaling
            Signedness: Signed
            WordLength: 50
        FractionLength: 35

Control Number of Jacobi Iterations

Use the numberOfSweeps parameter to control the number of iterations performed.

 m = 5;
 n = 3;
 rng('default');
 A = 10*randn(m,n)

A = 5×3

    5.3767  -13.0769  -13.4989
   18.3389   -4.3359   30.3492
  -22.5885    3.4262    7.2540
    8.6217   35.7840   -0.6305
    3.1877   27.6944    7.1474

Define fixed-point types that will never overflow. First, use the fixed.singularValueUpperBound function to determine the upper bound on the singular values. Then, define the integer length based on the value of the upper bound, with one additional bit for the sign, another additional bit for intermediate CORDIC growth, and one more bit for intermediate growth to compute the Jacobi rotations. Compute the fraction length based on the integer length and the desired word length.

svdUpperBound = fixed.singularValueUpperBound(m,n,max(abs(A(:))));
additionalBitGrowth = 3;
integerLength = ceil(log2(svdUpperBound)) + additionalBitGrowth;
wordLength = 16;
fractionLength = wordLength - integerLength;

Cast the matrix A to the signed fixed-point type.

T.A = fi([],1,wordLength,fractionLength,'DataType','Fixed');
A = cast(A,'like',T.A)

A = 
    5.3750  -13.0625  -13.5000
   18.3438   -4.3438   30.3438
  -22.5938    3.4375    7.2500
    8.6250   35.7812   -0.6250
    3.1875   27.6875    7.1562

          DataTypeMode: Fixed-point: binary point scaling
            Signedness: Signed
            WordLength: 16
        FractionLength: 5

Find the singular value decomposition of the fixed-point matrix A and calculate the relative error.

[U,s,V] = fixed.jacobiSVD(A);
relativeError = norm(double(U*diag(s)*V' - A))/norm(double(A))

relativeError = 
0.0069

Recalculate the singular value decomposition using only two Jacobi iterations and compute the relative error.

[U,s,V] = fixed.jacobiSVD(A,1);
relativeError = norm(double(U*diag(s)*V' - A))/norm(double(A))

relativeError = 
0.0704

Compare the relative error between the default 10 iterations and 2 iterations.

Input Arguments

A — Input matrix

matrix

Input matrix, specified as a matrix. A can be a signed fixed-point fi, a signed scaled double fi, double, or single data type.

Data Types: single | double | fi
Complex Number Support: Yes

numberOfSweeps — Number of Jacobi iterations

10 (default) | positive integer-valued scalar

Number of Jacobi iterations, specified as a positive integer-valued scalar. Increasing the number of Jacobi iterations improves the accuracy of the singular value decomposition.

Data Types: single | double | int8 | int16 | int32 | int64 | uint8 | uint16 | uint32 | uint64 | fi

Output Arguments

U — Left singular vectors

columns of matrix

Left singular vectors, returned as the columns of a matrix.

For fixed-point and scaled-double inputs, U is returned as a signed fixed-point or scaled-double fi with the same word length as A and fraction length equal to two less than the word length. One of these integer bits is used for the sign. The other integer bit allows +1 to be represented exactly.

For floating-point input, U has the same data type as A.

s — Singular values

column vector

Singular values, returned as a column vector. The singular values are nonnegative and returned in decreasing order. The singular values S have the same data type as A.

V — Right singular vectors

columns of matrix

Right singular vectors, returned as the columns of a matrix.

For fixed-point input and scaled-double inputs, V is returned as a signed fixed-point or scaled-double fi with the same word length as A and fraction length equal to two less than the word length. One of these integer bits is used for the sign. The other integer bit allows +1 to be represented exactly.

For floating-point input, V has the same data type as A. One of these integer bits is used for the sign, and the other integer bit is so that +1 can be represented exactly.

Limitations

If the input matrix A is not full rank and has more rows than columns, then the output U matrix is orthonormal up to the rank of A. $$A=U*diag(s)*V\text{'}$$ is still valid and V is orthonormal.

For example, if $$A=\left[\begin{array}{ccc}1& 1& 1\\ 2& 2& 2\\ 3& 3& 3\\ 0.5& 2& 4\end{array}\right]$$, then A has rank 2 and $$U\text{'}\ast U=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 0.8940\end{array}\right]$$.

If the input matrix A is not full rank and has more columns than rows, then the output V matrix is orthonormal up to the rank of A. $$A=U*diag(s)*V\text{'}$$ is still valid and U is orthonormal.

For example, if $$A\text{ =}\left[\begin{array}{cccc}1& 1& 1& 1\\ 2& 2& 2& 2\\ 1& 2& 3& 4\end{array}\right]$$, then A has rank 2 and $$V\text{'}\ast V=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 0.9653\end{array}\right]$$.

If the input matrix is full rank, then $$U\text{'}*U=I$$, where I = eye(n).

Tips

• The fixed.jacobiSVD function generates an economy sized vector output of the singular value decomposition. [U,s,V] = fixed.jacobiSVD(A) produces a vector s and unitary matrices U and V such that the dimensions of U, s and V are the same as the dimensions of svd with the "econ" and "vector" flags: [U,s,V] = svd(A,"econ","vector").

• The behaviors of the Square Jacobi SVD HDL Optimized and Non-Square Jacobi SVD HDL Optimized blocks are equivalent to [U,s,V] = fixed.jacobiSVD(A) when the input matrix A is square or non-square, respectively. If the input data type is fixed point with binary-point scaling, the function and the block provide bit-exact results. However, if the input data type is floating point, small numerical differences may exist between the function and the block.

  The numberOfSweeps input argument for the fixed.jacobiSVD function is equivalent to the Number of Jacobi iterations block parameter for both blocks.

Algorithms

The fixed.jacobiSVD function uses the two-sided Jacobi algorithm for singular value decomposition (SVD) [1][2][3]. Compared to the sequential Golub-Kahan-Reinsch algorithm for SVD [4], the Jacobi algorithm has inherent parallelism and performs better for FPGA and ASIC applications [5]. The Jacobi method is an iterative algorithm. The numberOfSweeps parameter determines the number of iterations performed. Most sources indicate that 10 iterations is sufficient for the Jacobi algorithm to converge.

References

[1] Jacobi, Carl G. J. “Über ein leichtes Verfahren die in der Theorie der Säcularstörungen vorkommenden Gleichungen numerisch aufzulösen.” Journal fur die reine und angewandte Mathematik 30 (1846): 51–94.

[2] Forsythe, George E., and Peter Henrici. “The Cyclic Jacobi Method for Computing the Principal Values of a Complex Matrix.” Transactions of the American Mathematical Society 94, no. 1 (January 1960): 1–23. https://doi.org/10.1090/S0002-9947-1960-0109825-2.

[3] Shiri, Aidin and Ghader Khosroshahi. 2019. “An FPGA Implementation of Singular Value Decomposition.” ICEE 2019:27th Iranian Conference on Electrical Engineering, Yazd, Iran, April 30–May 2, 2019, 416–22. IEEE. https://doi.org/10.1109/IranianCEE.2019.8786719.

[4] Golub, Gene H., and Charles F. Van Loan. Matrix Computations, 4th ed. Baltimore, MD: Johns Hopkins University Press, 2013.

[5] Athi, Mrudula V., Seyed R. Zekavat, and Allan A. Struthers. “Real-Time Signal Processing of Massive Sensor Arrays via a Parallel Fast Converging SVD Algorithm: Latency, Throughput, and Resource Analysis.” IEEE Sensors Journal 16, no. 8 (January 2016): 2519–26.https://doi.org/10.1109/JSEN.2016.2517040.

Extended Capabilities

C/C++ Code Generation

Generate C and C++ code using MATLAB® Coder™.

The fixed.jacobiSVD function generates efficient, purely integer C code.

Version History

Introduced in R2023a

R2024b: Bit-exact results with Jacobi singular value decomposition Simulink blocks

If the input data type is fixed point with binary-point scaling, fixed.JacobiSVD now provides bit-exact results with the Square Jacobi SVD HDL Optimized and Non-Square Jacobi SVD HDL Optimized blocks. Small numerical differences may exist compared to previous releases.