Main Content

fixed.jacobiSVD

Fixed-point Jacobi singular value decomposition

Since R2023a

Description

example

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

example

[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.

example

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

Examples

collapse all

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.4502
    1.6904
    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 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.1583    0.2715
   -0.6411   -0.7532
    0.7036   -0.5066
   -0.2615    0.3201

          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.1267
    0.1267    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.3578    3.1882
   18.3815  -13.0908
  -22.5636   -4.3423
    8.6213    3.4633

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

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.0062

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.0711

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

Input Arguments

collapse all

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

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

collapse all

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.

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.

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.

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 behavior of the Square Jacobi SVD HDL Optimized block is equivalent to [U,s,V] = fixed.jacobiSVD(A). The fixed.JacobiSVD function uses the same algorithm as the Square Jacobi SVD HDL Optimized block, with the same output data types. However, small numerical differences may exist in the least significant bit between the function and the block.

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

Version History

Introduced in R2023a