Updating Singular Value Decomposition for Rank One Matrix Perturbation
Ratnik Gandhi Amoli Rajgor
cs.LG math.NA









ratnik.gandhi@ahduni.edu.in


amoli.rajgor@iet.ahduni.edu.in


School of Engineering & Applied Science, Ahmedabad University, Ahmedabad-380009, India


An efficient Singular Value Decomposition (SVD) algorithm is an important tool for distributed and streaming computation in big data problems. It is observed that update of singular vectors of a rank-1 perturbed matrix is similar to a Cauchy matrix-vector product. With this observation, in this paper, we present an efficient method for
updating Singular Value Decomposition of rank-1 perturbed matrix in O(n^2  log(1/ϵ)) time. The method uses 
Fast Multipole Method (FMM) for updating singular vectors in O(n  log (1/ϵ)) time, where ϵ is the precision of computation. 



Updating SVD Rank-1 perturbation Cauchy matrixFast Multipole Method









§ INTRODUCTION

	SVD is a matrix decomposition technique having wide range of applications, such as image/video compression, real time recommendation system, text mining (Latent Semantic Indexing (LSI)), signal processing, pattern recognition, etc. Computation  of  SVD  is rather simpler in case of centralized system where entire data matrix is available at a single location. It is more complex when the matrix is distributed over a network of devices. Further, increase in  complexity is attributed to the real-time arrival of new data. Processing of data over distributed stream-oriented systems require efficient algorithms that generate results and update them in real-time. In streaming environment the data is continuously updated and thus the output of any operation performed over the data must be updated accordingly. In this paper, a special case of SVD updating algorithm is presented where the updates to the existing data are of rank-1. Organization of the paper is as follows: Section <ref> motivates the problem of updating SVD for a rank-1 perturbation and presents a characterization to look at the problem from matrix-vector product point of view. An existing algorithm for computing fast matrix-vector product using interpolation is discussed in Section <ref>. Section <ref> introduces 
Fast Multipole Method (FMM).
In Section <ref> we present an improved algorithm based on FMM for rank-1 SVD update that runs in O(n^2 log1/ϵ) time, where real ϵ >0 is a desired accuracy parameter. Experimental results of the presented algorithms are given in Section <ref>.
 
For completeness we have explained in details matrix factorization (<ref>), solution to Sylvester Matrix (<ref>), FAST algorithm for Cauchy matrix vector product (<ref>) and Fast Multipole Method (<ref>) at the end of this paper.

§ RELATED WORK

In a series of work Gu and others <cit.> present work on rank-one approximate SVD update. Apart from the low-rank SVD update, focus of their work is to discuss numerical computations and related accuracies in significant details. This leads to accurate computation of singular values, singular vectors and Cauchy matrix-vector product. Our work differs from this work as follows:  we use matrix factorization  that is explicitly based on solution of Sylvester equation <cit.> to reach Eq. (<ref>) that computes updated singular vectors (see Section <ref>). Subsequently, we reach an equation, Eq. (<ref>) (similar to equation (3.3) in <cit.>). Further, we show that with this new matrix decomposition we reach the computational complexity O(n^2 log1/ϵ) for updating rank-1 SVD.




§ SVD OF RANK-ONE PERTURBED MATRICES


	
	Let SVD of a m × n matrix A = UΣ V^⊤.[We consider the model and characterization as in <cit.>. A detailed factorization of this matrix is given in <ref>.]  Where, U ∈ℝ^m× m, Σ∈ℝ^m× n and V ∈ℝ^n× n, where, without loss of generality, we assume, m ≤ n. Let there be a rank-one update ab^⊤ to matrix A given by, 
	
	
    Â  =  A + ab^⊤

	
	and let ÛΣ̂V̂^⊤ denote the new(updated) SVD, where a ∈ℝ^m, b ∈ℝ^n,
	




	Thus,
	
	
    ÂÂ^⊤ = ÛΣ̂Σ̂^⊤Û^⊤  .

	An algorithm for updating SVD of a rank-1 perturbed matrix 
			is given in Bunch and Nielsen <cit.>. The algorithm updates singular values using characteristic polynomial and computes the updated singular vectors explicitly using the updated singular values.
	From (<ref>),
	
    ÂÂ^⊤    =     (U Σ V^⊤  +ab^⊤ )(U Σ V^⊤  + ab^⊤ )^⊤
        =     U Σ V^⊤ V Σ^⊤  U^⊤  + U Σ V^⊤ b_b̃a^⊤  + ab^⊤ V Σ^⊤ U^⊤_b̃^⊤ + ab^⊤ b_βa^⊤
    ÂÂ^⊤    =     UΣΣ^⊤ U^⊤  + b̃a^⊤  + ab̃^⊤  + β aa^⊤

	
		Where, b̃ = U Σ V^⊤ b,  b̃^⊤  = b^⊤ V Σ^⊤ U^⊤ and β = b^⊤ b.

	From (<ref>) it is clear that the problem of rank-1 update (<ref>) is modified to problem of three rank-1 updates (that is further converted to two rank-1 updates in (<ref>)).
	From (<ref>) and (<ref>) we get,
	
	
	
    ÛΣ̂Σ̂^⊤_D̂Û^⊤    =     UΣΣ^⊤_DU^⊤  + b̃a^⊤  + ab̃^⊤  + β aa^⊤
    ÛD̂Û^⊤    =    UDU^⊤  + ρ_1a_1a_1^⊤_ŨD̃Ũ^⊤ + ρ_2b_1b_1^⊤ .

	Refer Appendix <ref> Eq. (<ref>) for more details.
	Similar computation for right singular vectors is required.

			The following computation is to be done for each rank-1 update,  i.e., the procedure below will repeat four times, two times each for updating left and right singular vectors. From (<ref>) we have,
			
    ŨD̃Ũ^⊤    =     UDU^⊤  + ρ_1a_1a_1^⊤
        =     U(D + ρ_1a̅a̅^⊤ )_BU^⊤ ,

			where a̅ = U^⊤ a_1.
			From (<ref>) we have,
			
    B     :=     D + ρ_1a̅a̅^⊤
    		
    			B     =    C̃D̃C̃^⊤ (Schur-decomposition).

			From (<ref>) and (<ref>) we get,
			
			
    ŨD̃Ũ^⊤   =  U(C̃_ŨD̃C̃^⊤ )U^⊤_Ũ^⊤.

			After adding the rank-1 perturbation to UDU^⊤ in  (<ref>), the updated singular vector matrix is given by matrix-matrix product
			
    Ũ = U C̃.
 
	
		Stange <cit.>, extending
		the work of <cit.>, presents an efficient way of updating SVD by exploring the underlying structure of the matrix-matrix computations of (<ref>). 

		
			
	

 §.§ Updating Singular Values
 
	An approach for computing singular values is through eigenvalues. Given a matrix Ŝ = S + ρ uu^⊤, its eigenvalues d̃ can be computed in O(n^2) numerical operations by solving characteristic polynomial of (S + ρ uu^⊤ )x = d̃x, where S= diag(d_i). 
	Golub <cit.> has shown that the above characteristic polynomial has following form. 
	
	
    w(d̃) = 1 + ρ∑_i=1^nu_i^2/d_i - d̃.
 
	
	
	






	Note that in the equation above d̃ is an unknown and thus, though the polynomial function is structurally similar to Eq. (<ref>), we can not use FMM for solving it.
	
	Recall, in order to compute singular values D̂ of updated matrix Â we need to update D twice as there are two symmetric rank-1 updates, i.e., for 
    ŨD̃Ũ^⊤  = UDU^⊤  + ρ_1a_1a_1^⊤
 
	we will update B = D + ρ_1a̅a̅^⊤ Eq. (<ref>) and similarly for 
	
    ÛD̂Û^⊤  = ŨD̃Ũ^⊤  + ρ_2b_1b_1^⊤
 
	we will update B_1 = D̃ + ρ_2b̅b̅^⊤. 

	
	
	At times, while computing eigen-system, some of the eigenvalues and eigenvectors are known (This happens when there is some prior knowledge available about the eigenvalues or when the eigenvalues are approximated by using methods such as power iteration.). In such cases efficient SVD update method should focus on updating unknown eigen values and eigen vectors. Matrix Deflation is the process of eliminating known eigenvalue from the matrix.
	
	Bunch, Nielsen and Sorensen <cit.> extended Golub's work by bringing the notion of deflation. They presented a computationally efficient method for computing singular values by deflating the system for cases: (1) when some values of a̅ and b̅ are zero, (2) a̅ = 1 and b̅ = 1 and (3) B and B_1 has eigenvalues with repetition (multiplicity more than one). After deflation,
	the singular values of the updated matrix D̂  can be obtained by Eq. (<ref>).
	






	

 §.§ Updating Singular Vectors

	In order to update singular vectors the matrix-matrix product of (<ref>) is required. A naive method for matrix multiplication has complexity O(n^3). We exploit matrix factorization in Stange <cit.> that shows the structure of matrix C̃ to be Cauchy.

	
	From (<ref>) and (<ref>) we get,	

    C̃D̃C̃^⊤    =     D + ρ_1a̅a̅^⊤
    
    
    DC̃ - C̃D̃    =     - ρ_1a̅a̅^⊤C̃.

Equation (<ref>) is Sylvester equation with solution,

    (I_n ⊗ D + (-D̃)^⊤⊗ I_n)C̃ =  (- ρ_1a̅a̅^⊤C̃).


Simplifying L.H.S. of (<ref>) for C̃ we get,

    C̃ = [
    [ a̅_1     ;         ⋱;           a̅_n ]]
    [
    [ (μ_1 - λ_1)^-1              … (μ_n - λ_1)^-1;              ⋮                             ⋮; (μ_1 - λ_n)^-1              … (μ_n - λ_n)^-1 ]]
    [
    [ a̅^⊤ c_1         ;                 ⋱;                   a̅^⊤ c_n ]].


Where, C̃ = [[ c_1   … c_n ]] and from R.H.S. of (<ref>), c_i ∈ℝ^n are  columns of the form, 

    c_i = ρ_1 [[ a̅_1/μ_i - λ_1;              ⋮; a̅_n/μ_i - λ_n ]] a̅^⊤ c_i.





Simplifying (<ref>) further we get,


    c_i     = ρ_1 [[ a̅_1/μ_i - λ_1;              ⋮; a̅_n/μ_i - λ_n ]] [[ a̅_1 a̅_2 ]] [[ c_i1; c_i2 ]]
    
    c_i     = -ρ_1 (a̅_1c_i1 +  a̅_2c_i2) [[ a̅_1/λ_1 - μ_i;              ⋮; a̅_n/λ_n - μ_i ]].

Thus, denoting α̂ as -ρ_1 (a̅_1c_i1 +  a̅_2c_i2) we get,



    c_i = α̂[[ a̅_1/λ_1 - μ_i  …  a̅_n/λ_n - μ_i ]]^⊤ .

 
Placing (<ref>) in (<ref>) and we have,

    α̂[[ a̅_1/λ_1 - μ_i;              ⋮; a̅_n/λ_n - μ_i ]]     = α̂ρ_1 [[ a̅_1/μ_i - λ_1;              ⋮; a̅_n/μ_i - λ_n ]] ∑_k = 1^na̅_k^2/λ_k - μ_i.

Therefore for each element of c_i,

    α̂a̅_j/λ_j - μ_i    = α̂ρ_1  a̅_j/μ_i - λ_j ∑_k = 1^na̅_k^2/λ_k - μ_i
    
    
    
    
    0    = α̂(1 + ρ_1 ∑_k = 1^na̅_k^2/λ_k - μ_i)_=0

Equation inside the bracket of (<ref>) is the characteristic equation (<ref>) used for finding updated eigenvalues μ_i of D̃ by using eigenvalues of D i.e. λ. As μ_i are the zeros of this equation. Placing λ_k and μ_i will make the term in the bracket (<ref>) zero. Due to independence of the choice of scalar α̂, any value of α̂ can be used to scale the matrix C̅. In order to make the final matrix orthogonal, each column of C̅ is scaled by inverse of Euclidean norm of the respective column (<ref>).


From Eq. (<ref>) matrix notation for C̃ can be written as


    C̃ = [
    [ a̅_1     ;         ⋱;           a̅_n ]]
    [[   1/λ_1 - μ_1  ⋯  1/λ_1-μ_n;                          ⋮⋮; 1/λ_n - μ_1  ⋯  1/λ_n - μ_n ]]^C_C̅ = [[ c_1   … c_n ]][
    [ |c_1|      ;           ⋱;             |c_n| ]]^-1.

Where, c_i is the Euclidean norm of c_i.
From (<ref>) it is evident that the C matrix is similar to Cauchy matrix and C̃ is the scaled version of Cauchy matrix C. In order to update singular vectors we need to calculate matrix-matrix product as given in  (<ref>). From (<ref>) and (<ref>) we get,

    UC̃    =     U[
    [ a̅_1     ;         ⋱;           a̅_n ]]
    [[   1/λ_1 - μ_1  ⋯  1/λ_1-μ_n;                          ⋮⋮; 1/λ_n - μ_1  ⋯  1/λ_n - μ_n ]] 
    [
    [ |c_1|      ;           ⋱;             |c_n| ]]^-1




















    UC̃    =    [[ û_1   … û_n ]] C [
    [ |c_1|      ;           ⋱;             |c_n| ]]^-1.
         Where, û_i     =     u_ia̅_i  and  i = 1, … ,n.
    
    UC̃    =     U_1  C [
    [ |c_1|      ;           ⋱;             |c_n| ]]^-1
    
    UC̃    =     U_2 [
    [ |c_1|      ;           ⋱;             |c_n| ]]^-1
         Where,  U_1     =    [[ û_1   … û_n ]] and  U_2 = U_1 C .
 

	

  §.§.§ Trummer's Problem: Cauchy Matrix-Vector Product

	
	In (<ref>) there are n vectors in U which are multiplied with Cauchy matrix C. The problem of multiplying a Cauchy matrix with a vector is called Trummer's problem. As there are n vectors in (<ref>), matrix-matrix product in it can be represented as n matrix-vector products, i.e., it is same as solving Trummer's problem n times. Section <ref> describes an algorithm given in <cit.> which efficiently computes such matrix-vector product in O(n  log^2(n)) time.
	  
	
	

§ FAST: METHOD BASED ON POLYNOMIAL INTERPOLATION AND FFT
 





























	Consider the the matrix - (Cauchy) matrix product U_2 = U_1C of Section <ref>. This product can be written as
	
	
	
	
    U_2     =    [[ û_1   … û_n ]] C
    
    	U_2     =    [[ u_11  ⋯   u_1n;            ⋮⋮⋮; u_n1  ⋯   u_nn ]][[   1/λ_1 - μ_1  ⋯  1/λ_1-μ_n;                         ⋮⋮⋮; 1/λ_n - μ_1  ⋯  1/λ_n - μ_n ]].

	The dot product of each row vector of U_1 and the coulmn vector of the Cauchy matrix can be represented in terms of a function of eigenvalues μ_i,  i = 1,…,n,
	
	
    f(μ_1)     =    [[ u_11  ⋯   u_1n ]][[ 1/λ_1 - μ_1;           ⋮; 1/λ_n - μ_1 ]]
        =    u_11/λ_1 -μ_1+…+u_1n/λ_n -μ_1
     
    	f(μ_1)    =    ∑_j=1^nu_1j/λ_j - μ_1.
    Hence in general,      
    
    	f(x)     =    ∑_j=1^nu_j/λ_j - x.

	Equation (<ref>) can be shown as the ratio of two polynomials 
	
    f(x)     =    h(x)/g(x),     where,   g(x) = ∏_j=1^n (λ_j - x).
    
    		h(x)    =    g(x)∑_j=1^nu_j/λ_j - x
    
    		h(x)    =   ∏_j=1^n (λ_j - x)∑_j=1^nu_j/λ_j - x


		The FAST algorithm  <cit.> represents the a matrix-vector product as (<ref>) and (<ref>). It then finds solutions to this problem by the use of interpolation.		
	The FAST algorithm of Gerasoulis <cit.> has been reproduce in <ref>.

	
	

§ MATRIX-VECTOR PRODUCT USING FMM
 
	The FAST algorithm  in <cit.> has complexity O (n  log^2n). It computes the function f(x) using FFT and interpolation. We observe that these two methods are two major procedures that contributes to the overall complexity of FAST algorithm. To reduce this complexity we present an algorithm that uses FMM for finding Cauchy matrix-vector product that updates the SVD of rank-1 perturbed matrix with time complexity  O(n^2  log(1/ϵ)).
	
	Recall from Section <ref>, update of singular vectors require matrix-(Cauchy) matrix product U_2 = U_1C, This product is represented as the function below.
	
    f(μ_i)  =  ∑_j=1^n-u_j/μ_i - λ_j.

	We use FMM to compute this function. 

	




	

 §.§ Fast Multipole Method (FMM)

	An algorithm for computing potential and force of n particles in a system is given in Greengard and Rokhlin <cit.>. This algorithm enables fast computation of forces in an n-body problem by computing interactions among points in terms of far-field and local expansions. It starts with clustering particles in groups such that the inter-cluster and intra-cluster distance between points make them well-separated. Forces between the clusters are computed using Multipole expansion. 

			
	Dutt et al. <cit.> describes the idea of FMM for particles in one dimension and presents an algorithm for evaluating sums of the form, 
						
	
    f(x) = ∑_k=1^Nα_k ·ϕ(x - x_k).

						
	Where, {α_1,…,α_N} is a set of complex numbers and ϕ(x) can be a function that is singular at x=0 and smooth every where. Based on the choice of function, (<ref>) can be used to evaluate sum of different forms for example:
										
	
    f(x) = ∑_k=1^Nα_k/x-x_k

	where, ϕ(x) = 1/x.
	Dutt et al. <cit.> presents an algorithm for computing  summation (<ref>) using FMM that runs in O(n  log(1/ϵ)). 
			








			 
			


 §.§ Summation using FMM

			
Consider well separated points {x_1,x_2,…, x_N} and {y_1,…,y_M} in ℝ such that for points x_0,y_0 ∈ℝ and r>0, r ∈ℝ.

    |x_i - x_0|     <     r  ∀ i = 1,…,N
    
    |y_i - y_0|     <     r  ∀ i = 1,…,M  and
    
    |x_0 - y_0|     >     4r

For a function defined as f:ℝ→ℂ such that
				
    f(x) = ∑_k=1^Nα_k/x-x_k

where, {α_1 , …, α_k} is a set of complex numbers. Given f(x), the task is to find f(y_1),…,f(y_M). 




§ RANK-ONE SVD UPDATE
	
In this section, we present Algorithm <ref> that uses FMM and updates Singular Value Decomposition in O(n^2  log(1/ϵ)) time. 
 



 §.§ Algorithm: Update SVD of rank-1 modified matrix using FMM



  INPUT A ∈ℝ^m × n = UΣ V^T, a∈ℝ^mand b∈ℝ^n


  OUTPUT U_n, Σ_n, and V_n


  STEP 1 Compute b̃ = UΣ V^Tb, ã = VΣ^TU^Ta, β = b^Tb, α = a^Ta, D_u = ΣΣ^T and D_v = Σ^TΣ.


  STEP 2 Compute the Schur decomposition of [[  β   1; 1    0 ]] = Q_u [[ ρ_1    0;  0   ρ_2 ]]Q_u^T and [[ a  b̃ ]]Q_u = [[ a_1  b_1 ]].

  STEP 3 Compute the Schur decomposition of [[  α   1; 1    0 ]] = Q_v [[ ρ_1    0;  0   ρ_2 ]]Q_v^T and [[ b  ã ]]Q_v = [[ a_2  b_2 ]].


  STEP 4 Compute updated left singular vector Ũ and D̃ by calling the procedure RankOneUpdate(U,a_1,D,ρ_1) - Algorithm <ref>.


  STEP 5 Compute left singular vector of Â, i.e., U_n  and D_n by calling the procedure RankOneUpdate(Ũ,b_1D̃,ρ_2).


  STEP 6 Compute updated right singular vector Ṽ and D̃_̃ṽ by calling the procedure RankOneUpdate(V,a_2,D_v,ρ_3).


  STEP 7 Compute right singular vector of Â i.e. V_n   and D_vn by calling the procedure RankOneUpdate(Ṽ,b_2,D̃_̃ṽ,ρ_4).


  STEP 8 Find singular values by computing square root of the updated eigenvalues Σ_n = √(D_n).

Note that the Schur decomposition of Steps 2 and 3 are computed over constant size matrices and thus the decomposition will take constant time.


 §.§ Algorithm: RankOneUpdate




  INPUT U,a_1,D,ρ_1

  OUTPUT Ũ, D̃


  STEP 1 Compute a̅ = U^T a_1. 






  STEP 2 Compute μ as zeros of the equation w(μ) = 1 + ρ_1 ∑_i=1^na̅_i^2/λ_i - μ.


  STEP 3 Compute C = [[   1/λ_1 - μ_1  ⋯  1/λ_1-μ_n;                          ⋮⋮; 1/λ_n - μ_1  ⋯  1/λ_n - μ_n ]].


  STEP 4 Compute U_1 = U [[ a̅_1     ;    ⋱;           a̅_n ]]



  STEP 5Compute C̅ = [[ a̅_1     ;    ⋱;           a̅_n ]] [[   1/λ_1 - μ_1  ⋯  1/λ_1-μ_n;                          ⋮⋮; 1/λ_n - μ_1  ⋯  1/λ_n - μ_n ]]


  STEP 6 Compute U_2= U_1C as n matrix-vector product. Where each row-column dot product is represented as a function 

    f(μ_i)  =  ∑_j=1^n-u_j/μ_i - λ_j .
 for each μ_i.

U_2 = FMM(λ, μ, -U_1) - (<ref>).


  STEP 7 Form Ũ by dividing each column of U_2 by norm of respective column of C̅.



	





































	


 §.§ Complexity for Rank-One SVD update

	theoremTheorem
			
			Given a matrix A ∈ℝ^m × n such that m ≤ n and precision of computation parameter ϵ>0, the complexity of computing SVD of a rank-1 update of A with Algorithm <ref> is O(n^2  log (1/ϵ)).
			
			
			
			Follows from Algorithm <ref> and Table <ref>.
			
					
		
		
		





§ NUMERICAL RESULTS

All the computations for the algorithm were performed on MATLAB over a machine with Intel i5, quad-core, 1.7 GHz, 64-bit processor with 8 GB RAM. Matrices used in these experiments are square and generated randomly with values ranging from [1,9]. The sample size varies from 2 × 2 to 35 × 35. 
For computing the matrix vector product we use FMM Algorithm (instead of FAST Algorithm) with machine precision ϵ  =5^-10. 
















In order to update singular vectors we need to perform two rank-1 updates (<ref>). For each such rank-1 update we use FMM. 


    ÛD̂Û^T     =  UDU^T + ρ_1a_1a_1^T_ŨD̃Ũ^T + ρ_2b_1b_1^T.

Figure <ref> shows the time taken by FAST Algorithm and FMM Algorithm to compute the first rank-1 update. 







 §.§ Choice of ϵ

In earlier computations we fixed the machine precision parameter ϵ for FMM based on the order of the Chebyshev polynomials, i.e., ϵ = 5^-p, where p is the order of Chebyshev polynomial. As we increase the order of polynomials - we expect reduction in the error and increase in updated vectors accuracy. This increase in accuracy comes with the cost of higher computational time. To show this we fix the input matrix dimension to 25 × 25 and generate values of these matrices randomly in the range [0,1]. Figure <ref> shows the error between updated singular vectors generated by Algorithm <ref> and exact computation of singular vector updates. 
The results in figure justifies our choice of fixed machine precision p = 20.



Error is computed using the equation (<ref>) <cit.>, where Â is the perturbed matrix, ÛΣ̂V̂^T is the approximation computed using FMM-SVDU and max σ̂ is the directly computed maximum eigenvalue of Â. 


    Error    = max|Â - ÛΣ̂V̂^T/max σ̂|


Table <ref> shows the accuracy of rank-1 SVD update using FMM- SVDU for varying sample size. Figure <ref> shows plot of the accuracy of FMM-SVDU with varying sample size. 

	
		

		








		

§ CONCLUSION

		In this paper we considered the problem of updating Singular Value Decomposition of a rank-1 perturbed matrix. We presented an efficient algorithm for updating SVD of rank-1 perturbed matrix that uses the Fast Multipole method for improving Cauchy matrix- vector product computational efficiency. 

An interesting and natural extension of this work is to consider updates of rank-k.



				






§ REFERENCES




10
url<#>1urlprefixURL href#1#2#2 #1#1

gu93
M. Gu, Studies in numerical linear algebra, Ph.D. thesis, Yale University
  (1993).

ge93
M. Gu, S. C. Eisenstat, A stable and fast algorithm for updating the singular
  value decomposition, Tech. Rep. YALEU/DCS/RR-966, Dept. of Computer Science,
  Yale University (1993).

ge94
M. Gu, S. C. Eisenstat, A stable and efficient algorithm for the rank-one
  modification of the symmetric eigenproblem, SIAM journal on Matrix Analysis
  and Applications 15 (4) (1994) 1266–1276.

ge95
M. Gu, S. C. Eisenstat, Downdating the singular value decomposition, SIAM
  Journal on Matrix Analysis and Applications 16 (3) (1995) 793–810.

stange08
P. Stange, On the efficient update of the singular value decomposition, PAMM
  8 (1) (2008) 10827–10828.

bn78
J. R. Bunch, C. P. Nielsen, Updating the singular value decomposition,
  Numerische Mathematik 31 (2) (1978) 111–129.

golub73
G. H. Golub, Some modified matrix eigenvalue problems, Siam Review 15 (2)
  (1973) 318–334.

bnc78
J. R. Bunch, C. P. Nielsen, D. C. Sorensen, Rank-one modification of the
  symmetric eigenproblem, Numerische Mathematik 31 (1) (1978) 31–48.

gerasoulis88
A. Gerasoulis, A fast algorithm for the multiplication of generalized Hilbert
  matrices with vectors, Mathematics of Computation 50 (181) (1988) 179–188.

gr87
L. Greengard, V. Rokhlin, A fast algorithm for particle simulations, Journal of
  computational physics 73 (2) (1987) 325–348.

dgr96
A. Dutt, M. Gu, V. Rokhlin, Fast algorithms for polynomial interpolation,
  integration, and differentiation, SIAM Journal on Numerical Analysis 33 (5)
  (1996) 1689–1711.









§ MATRIX FACTORIZATION

We consider the matrix factorization in <cit.>. Let SVD of a m × n matrix A = UΣ V^⊤. Where, U ∈ℝ^m× m, Σ∈ℝ^m× n and V ∈ℝ^n× n. Here it is assumed that m ≤ n. Let there be a rank-one update ab^⊤ to matrix A Eq. (<ref>) and let ÛΣ̂V̂^⊤ denote the new(updated) SVD, where a ∈ℝ^m, b ∈ℝ^n.
	
	
    Â    =     A + ab^⊤
    Â    =    ÛΣ̂V̂^⊤

	
	In order to find SVD of this updated matrix Â we need to compute ÂÂ^⊤ and Â^⊤Â because left singular vector Û of Â is orthonormal eigenvector of ÂÂ^⊤ and right singular vector V̂ of Â is orthonormal eigenvector of Â^⊤Â.
	
	
    ÂÂ^⊤    =     (ÛΣ̂V̂^⊤ )(ÛΣ̂V̂^⊤ )^⊤
        =    ÛΣ̂V̂^⊤V̂Σ̂^⊤Û^⊤
    ÂÂ^⊤    =    ÛΣ̂Σ̂^⊤Û^⊤     Where, V̂^⊤V̂ = I

	
	
    ÂÂ^⊤    =     (U Σ V^⊤  +ab^⊤ )(U Σ V^⊤  + ab^⊤ )^⊤
        =     (U Σ V^⊤  +ab^⊤ )(V Σ^⊤  U^⊤  + ba^⊤ )
        =     U Σ V^⊤ V Σ^⊤  U^⊤  + U Σ V^⊤ b_b̃a^⊤  + ab^⊤ V Σ^⊤ U^⊤_b̃^⊤ + ab^⊤ b_βa^⊤
    ÂÂ^⊤    =     UΣΣ^⊤ U^⊤  + b̃a^⊤  + ab̃^⊤  + β aa^⊤     Where,  V^⊤ V = I

	
		Where, b̃ = U Σ V^⊤ b,  b̃^⊤  = b^⊤ V Σ^⊤ U^⊤ and β = b^⊤ b.

	From Eq. (<ref>) it is clear that the problem of rank-1 update Eq. (<ref>) is modified to problem of three rank-1 updates that is further converted to two rank-1 updates in Eq. (<ref>).
	Equating Eq. (<ref>) and Eq. (<ref>) we get,
	
	
	
    ÛΣ̂Σ̂^⊤_D̂Û^⊤    =     UΣΣ^⊤_DU^⊤  + b̃a^⊤  + ab̃^⊤  + β aa^⊤
    ÛD̂Û^⊤    =     UDU^⊤ + b̃a^⊤  + ab̃^⊤  + β aa^⊤
        =     UDU^⊤ + β aa^⊤  + b̃a^⊤  + ab̃^⊤
        =     UDU^⊤ + (β a + b̃)a^⊤  + ab̃^⊤
        =     UDU^⊤ + [[ β a + b̃   a ]]  [[  a^⊤; b̃^⊤ ]]

	
    ÛD̂Û^⊤    =     UDU^⊤  + [[ a   b̃ ]] [[  β   1; 1    0 ]] [[  a^⊤; b̃^⊤ ]]

	
    ÛD̂Û^⊤    =     UDU^⊤  + [[ a   b̃ ]] Q [[ ρ_1    0;  0   ρ_2 ]]Q^⊤[[  a^⊤; b̃^⊤ ]]

	
    Where,       [[  β   1; 1    0 ]]= Q [[ ρ_1    0;  0   ρ_2 ]]Q^⊤ (Schur-decomposition)

	
    ÛD̂Û^⊤    =     UDU^⊤  + [[ a   b̃ ]] Q_[[ a_1    b_1 ]][[ ρ_1    0;  0   ρ_2 ]]Q^⊤[[  a^⊤; b̃^⊤ ]]_[[ a_1    b_1 ]]^⊤
    ÛD̂Û^⊤    =     UDU^⊤  + [[ a_1    b_1 ]] [[ ρ_1    0;  0   ρ_2 ]][[ a_1^⊤; b_1^⊤ ]]^⊤
    ÛD̂Û^⊤    =    UDU^⊤  + ρ_1a_1a_1^⊤_ŨD̃Ũ^⊤ + ρ_2b_1b_1^⊤ .

		Similarly for computing right singular vectors do the following.
		
    Â^⊤Â    =     (U Σ V^⊤  +ab^⊤ )^⊤ (U Σ V^⊤  + ab^⊤ ) 
        =     (V Σ^⊤  U^⊤  + ba^⊤ )(U Σ V^⊤  +ab^⊤ )
    Â^⊤Â    =     V Σ^⊤  U^⊤ U Σ V^⊤  + V Σ^⊤  U^⊤ a_ãb^⊤  + ba^⊤ U Σ V^⊤_ã^⊤+ ba^⊤ a_αb^⊤
    
    			(V̂Σ̂^⊤Û^⊤ )(ÛΣ̂V̂^⊤ )    =     VΣ^⊤Σ V^⊤  + ãb^⊤  + bã^⊤  + α bb^⊤     Where,  U^⊤ U = I .
    V̂Σ̂^⊤Σ̂V̂^⊤    =      VΣ^⊤Σ V^⊤  + ρ_3a_2a_2^⊤  + ρ_4b_2b_2^⊤

			
			The following computation is to be done for each rank-1 update in Eq. (<ref>) and Eq. (<ref>) i.e. the below procedure will repeat four times, two times each for updating left and right singular vectors. From Eq. (<ref>) we have,
			
    ŨD̃Ũ^⊤    =     UDU^⊤  + ρ_1a_1a_1^⊤
        =     UDU^⊤  + ρ_1(UU^⊤ a_1)(UU^⊤ a_1)^⊤
        =     UDU^⊤  + ρ_1(Ua̅)(Ua̅)^⊤ Where, a̅ = U^⊤ a_1
        =     UDU^⊤  + ρ_1(Ua̅a̅^⊤ U^⊤ )
    ŨD̃Ũ^⊤    =     U(D + ρ_1a̅a̅^⊤ )_BU^⊤ .

			From Eq. (<ref>) we have,
			
    B     :=     D + ρ_1a̅a̅^⊤
    		
    			B     =    C̃D̃C̃^⊤ (Schur-decomposition).

			Placing Eq. (<ref>) in Eq. (<ref>) we get,
			
			
    ŨD̃Ũ^⊤   =  U(C̃_ŨD̃C̃^⊤ )U^⊤_Ũ^⊤.

			After the rank-1 update to UDU^⊤ Eq. (<ref>) the updated singular vector matrix is given by matrix-matrix product
			
    Ũ = U C̃.
 





§ SOLUTION TO SYLVESTER EQUATION

In Section <ref> we discussed method for updating singular vectors. For the same we derived solutions to (<ref>) using Sylvester equation. In this section we present details of how this solution can be obtained.


For simplicity consider the case where the dimension of all the matrices (C, D and D̃) is 2 × 2.


    L.H.S     = (I_n ⊗ D + (-D̃)^⊤⊗ I_n)C̃
        = {([[ 1 0; 0 1 ]] ⊗[[ λ_1   0;   0 λ_2 ]])  - ([[ μ_1   0;   0 μ_2 ]]^⊤⊗[[ 1 0; 0 1 ]])}[[ c_11 c_12; c_21 c_22 ]] 
        = {[[ λ_1   0   0   0;   0 λ_2   0   0;   0   0 λ_1   0;   0   0   0 λ_2 ]] - [[ μ_1   0   0   0;   0 μ_1   0   0;   0   0 μ_2   0;   0   0   0 μ_2 ]]}[[ c_11; c_21; c_12; c_22 ]] 
        = [[ λ_1 - μ_1         0         0         0;         0 λ_2 - μ_1         0         0;         0         0 λ_1 - μ_2         0;         0         0         0 λ_2 - μ_2 ]] [[ c_11; c_21; c_12; c_22 ]]

 
 
    R.H.S     =  (- ρ_1a̅a̅^⊤C̃)
        = -ρ_1(  [[ a̅_1; a̅_2;  ]][[ a̅_1 a̅_2 ]][[ c_11 c_12; c_21 c_22 ]])
        = -ρ_1(  [[ a̅_1; a̅_2;  ]][[ a̅_1c_11 + a̅_2 c_21 a̅_1c_12 + a̅_2 c_22 ]])
        = -ρ_1(  [[ a̅_1; a̅_2;  ]][[[ a̅_1 a̅_2 ]][[ c_11; c_21 ]]   [[ a̅_1 a̅_2 ]][[ c_12; c_22 ]]])
        = -ρ_1(  [[ a̅_1; a̅_2;  ]][[ a̅^⊤ c_1 a̅^⊤ c_2 ]])
       Where, C̃  = [[ c_1 c_2 ]] and   c_1 = [[ c_11; c_21 ]] c_2 = [[ c_12; c_22 ]]
        = -ρ_1(  [[ a̅_1a̅^⊤ c_1 a̅_1a̅^⊤ c_2; a̅_2a̅^⊤ c_1 a̅_2a̅^⊤ c_2 ]])
        = -ρ_1[[ a̅_1a̅^⊤ c_1; a̅_2a̅^⊤ c_1; a̅_1a̅^⊤ c_2; a̅_2a̅^⊤ c_2 ]]


Equating L.H.S and R.H.S we get,


    [[ λ_1 - μ_1         0         0         0;         0 λ_2 - μ_1         0         0;         0         0 λ_1 - μ_2         0;         0         0         0 λ_2 - μ_2 ]] [[ c_11; c_21; c_12; c_22 ]]     = -ρ_1[[ a̅_1a̅^⊤ c_1; a̅_2a̅^⊤ c_1; a̅_1a̅^⊤ c_2; a̅_2a̅^⊤ c_2 ]].

Therefore,

    [[ c_11; c_21; c_12; c_22 ]]     = -ρ_1[[ λ_1 - μ_1         0         0         0;         0 λ_2 - μ_1         0         0;         0         0 λ_1 - μ_2         0;         0         0         0 λ_2 - μ_2 ]]^-1[[ a̅_1a̅^⊤ c_1; a̅_2a̅^⊤ c_1; a̅_1a̅^⊤ c_2; a̅_2a̅^⊤ c_2 ]]
    [[ c_11; c_21; c_12; c_22 ]]     = -ρ_1[[ 1/λ_1 - μ_1           0           0           0;           0 1/λ_2 - μ_1           0           0;           0           0 1/λ_1 - μ_2           0;           0           0           0 1/λ_2 - μ_2 ]][[ a̅_1a̅^⊤ c_1; a̅_2a̅^⊤ c_1; a̅_1a̅^⊤ c_2; a̅_2a̅^⊤ c_2 ]]
    [[ c_11; c_21; c_12; c_22 ]]     = -ρ_1[[ a̅_1a̅^⊤ c_1/λ_1 - μ_1; a̅_2a̅^⊤ c_1/λ_2 - μ_1; a̅_1a̅^⊤ c_2/λ_1 - μ_2; a̅_2a̅^⊤ c_2/λ_2 - μ_2 ]]
    [[ c_11; c_21; c_12; c_22 ]]     = ρ_1[[ a̅_1a̅^⊤ c_1/μ_1 - λ_1; a̅_2a̅^⊤ c_1/μ_1 - λ_2; a̅_1a̅^⊤ c_2/μ_2 - λ_1; a̅_2a̅^⊤ c_2/μ_2 - λ_2 ]]
    Where, [[ c_11; c_21 ]]     = ρ_1[[ a̅_1a̅^⊤ c_1/μ_1 - λ_1; a̅_2a̅^⊤ c_1/μ_1 - λ_2 ]]   and [[ c_12; c_22 ]]  = ρ_1[[ a̅_1a̅^⊤ c_2/μ_2 - λ_1; a̅_2a̅^⊤ c_2/μ_2 - λ_2 ]]
    i.e.  c_1     = ρ_1[[ a̅_1/μ_1 - λ_1; a̅_2/μ_1 - λ_2 ]]a̅^⊤ c_1   and  c_2 = ρ_1[[ a̅_1/μ_2 - λ_1; a̅_2/μ_2 - λ_2 ]]a̅^⊤ c_2 
    In general c_i     = ρ_1[[ a̅_1/μ_i - λ_1;              ⋮; a̅_2/μ_i - λ_n ]]a̅^⊤ c_1


By placing c_1 and c_2 in C̃ = [[ c_1 c_2 ]] we get C̃ as below.


    C̃    = [[ c_1 c_2 ]]
    C̃    = [[ a̅_1a̅^⊤ c_1/μ_1 - λ_1 a̅_1a̅^⊤ c_2/μ_2 - λ_1; a̅_1a̅^⊤ c_1/μ_1 - λ_2 a̅_1a̅^⊤ c_2/μ_2 - λ_2 ]]
    C̃    = [[ a̅_1    0;    0 a̅_2 ]][[ a̅^⊤ c_1/μ_1 - λ_1 a̅^⊤ c_2/μ_2 - λ_1; a̅^⊤ c_1/μ_1 - λ_2 a̅^⊤ c_2/μ_2 - λ_2 ]]
    C̃    = [[ a̅_1    0;    0 a̅_2 ]][[ 1/μ_1 - λ_1 1/μ_2 - λ_1; 1/μ_1 - λ_2 1/μ_2 - λ_2 ]][[ a̅^⊤ c_1        0;        0 a̅^⊤ c_2 ]]



			



§ FAST ALGORITHM

	
In this section we present FAST Algorithm <cit.> that computes functions of the form (<ref>) using polynomial interpolation in time O(n log^2 n).
	
	
  * Compute the coefficients of g(x) in its power form, by using FFT polynomial multiplication, in O(n(log n)^2) time.


	
	Decription: Uses FFT for speedy multiplication which reduces complexity of multiplication to O(n  log(n)) from O(n^2).

	Input: Function g(x) and eigenvalues of D̃ i.e. [μ_1,…,μ_n]

	Output: Coefficients of g(x) i.e. [a_0,a_1,…,a_n]

	Complexity: O(n  log^2n)
	
	
  * Compute the coefficients of g'(x) in O(n) time.


	
	Description: Differentiate the function g(x) and then compute its coefficients.

	Input: Function g(x)

	Output: Coefficients of g'(x) i.e. [b_0,b_1,…,b_n]

	Complexity: O(n)

	
	
  * Evaluate g(λ_i), g'(λ_i) and g(μ_i).


	
	Description: For the functions g(x) and g'(x) evaluate their values at λ_i and μ_i

		Input: Function g(x) and g'(x), eigenvalues [λ_1,λ_2,…,λ_n] of D and [μ_1,μ_2,…,μ_n] of D̃

		Output: g(λ_i), g'(λ_i) and g(μ_i)

		Complexity: O(n  log^2n)

	
	
  * Compute h_j = u_jg'(λ_j).


	
	Description: For each eigenvalue λ_j compute its function value h(λ_j) i.e. find the points (λ_j,u_jg'(λ_j))

		Input: [λ_1,λ_2,…,λ_n]

		Output: h_j

		Complexity: O(n)

	
	
  * Find interpolation polynomial h(x) for the points (λ_j,h_j).


	
	Description: Given the function values and input i.e. points (λ_j,h_j) find interpolation polynomial for n points.

		Input: Points (λ_j,h_j)

		Output: h(x) 

		Complexity: O(n  log^2n)

	
	
  * Compute v_i = h_(μ_i)/g(μ_i).


	
	Description: Compute the ratio v_i = h_(μ_i)/g(μ_i) for each μ_i where, v_i is value of f(μ_i) at each μ_i.

		Input: Function h(x) and g(x) and [μ_1,μ_2,…,μ_n]

		Output:f(μ_i) = v_i

		Complexity: O(n)

	
	





§ FAST MULTIPOLE METHOD



 §.§ Interpolation and Chebyshev Nodes


				
Chebyshev nodes are the roots of the Chebyshev polynomials and they are used as points for interpolation. These nodes lie in the range [-1,1]. A polynomial of degree less than or equal to n-1 can fit over n Chebyshev nodes. For Chebyshev nodes the approximating polynomial is computed using Lagrange interpolation. 

Expansions are used to quantify interactions among points. Expansions are only computed for points which are well-separated from each other.



				

 §.§ FMM Algorithm

					
						
  STEP 1Description: Decide the size of Chebyshev expansion i.e., p = -log_5(ϵ) = log_5 (1/ϵ) where ϵ > 1 is the precision of computation or machine accuracy parameter.

						Input: ϵ

						Output: p

						Complexity: O(1)

						
						
  STEP 2Description: Set s as the number of points in the cell of finest level (s ≈ 2p) and level of finest division nlevs = log_2(N/s) where, N is the number of points.
		
						Input: N and p

						Output: s and nlevs

						Complexity: O(1)
	
						
						
  STEP 3Description: Consider p Chebyshev nodes defined over an interval [-1,1] of the form below.
						
						
    t_i = cos ( 2i - 1/p·π/2)

						Where, i = 1,…,p.

					
												
						
  STEP 4Description: Consider Chebyshev polynomials of the form 
						

    u_j(t) = ∏_k=1 
     k≠ j^pt - t_k/t_j - t_k.
 
Where, j = 1,…,p.

		

	
  STEP 5 Description: Calculate the far-field expansion Φ_l,i at i^th subinterval of level l.
						
	For i^th interval of level l (nlevs) far-field expansion due to points in interval [x_0 - r, x_0 + r] about center of the interval x_0 is defined by a vector of size p as 
						
    Φ_nlevs,i = ∑_k=1^Nα_k ·t_i/3r - t_i(x_k -x_0).

						Input: l = nlevs, α_k, x_k, x_0, r, t_i for i = 1,…,2^nlevs, 

						Output: Φ_nlevs,i

						Complexity: O(Np)

						
							
  STEP 6[Bottom-up approach]
 Description: Compute the far-field expansion of individual subintervals in terms of far-field expansion of their children. These are represented by p × p matrix defined as below.
												
							
    M_L(i,j)     =     u_j(t_i/2+t_i)
    
    							M_R(i,j)     =     u_j(t_i/2-t_i)

				Where, i = 1,…,p and j = 1,…,p. {u_1,…,u_p} and {t_1,…,t_p} are as defined in Eq. (<ref>) and Eq. (<ref>) respectively.
				
				Far-field expansion for i^th subinterval of level l due to far-field expansion of its children is computed by shifting children's far-field expansion by M_L or M_R and adding those shifts as below. 
    Φ_l,i = M_L ·Φ_l+1,2i-l + M_R ·Φ_l+1,2i

				Input: M_L, M_R and Φ

				Output: Φ_l,i

				Complexity: O(2Np^2/s)

				
				
  STEP 7 Description: Calculate the local expansion Ψ_l,i at i^th subinterval of level l.
									
				For i^th interval of level l local expansion due to points outside the interval [y_0 - r, y_0 + r] about center of the interval y_0 is defined by a vector of size p as
				
    Ψ_i = ∑_k=1^Nα_k ·1/rt_i - (x_k - x_0).

				Input: t_i, α_k, r, x_k,and  x_0

				Output: Ψ_i

				Complexity: O(Np)

				
				
  STEP 8[Top-down approach]
 Description: Compute the local expansion of individual subintervals in terms of local expansion of their parents. These are represented by p × p matrix defined as below. 
						
				
    S_L(i,j)     =     u_j(t_i - 1/2)
    
    				S_R(i,j)     =     u_j(t_i + 1/2)

				Where, i = 1,…,p and j = 1,…,p. {u_1,…,u_p} and {t_1,…,t_p} are as defined in Eq. (<ref>) and Eq. (<ref>) respectively.
				
				Compute local expansions using far-field expansion using p × p matrix defined as below.
						
				
    T_1(i,j)     =     u_j(3/t_i - 6)
    
    				T_2(i,j)     =     u_j(3/t_i - 4)
    
    				T_3(i,j)     =     u_j(3/t_i + 4)
    
    				T_4(i,j)     =     u_j(3/t_i + 6)

				
				Local expansion for each subinterval at finer level is computed using local expansion of their parents. For this, first the local expansion of parent interval are shifted by S_L or S_R and then the result is added with interactions of subinterval with other well separated subintervals (which were not considered at the parent level).
				
				
    Ψ_l+1,2i-1   =  S_L ·Ψ_l,i + T_1 ·Φ_l+1,2i-3 + T_3 ·Φ_l+1,2i+1 + T_4 ·Φ_l+1,2i+2
	
				
				
    Ψ_l+1,2i   =  S_R ·Ψ_l,i + T_1 ·Φ_l+1,2i-3 + T_2 ·Φ_l+1,2i-2 + T_4 ·Φ_l+1,2i+2

				Input: S_L, S_R, T_1, T_2, T_3, T_4, Φ and Ψ

				Output: Ψ_l+1,2i-1 and Ψ_l+1,2i

				Complexity: O(8Np^2/s)

				
				
  STEP 9 Description: Evaluate local expansion Ψ_nlevs,i at some of {y_j} (which falls in subinterval i of level nlevs) to obtain a p size vector. 

				Complexity: O(Np)

				
				
  STEP 10 Description: Add all the remaining interactions which are not covered by expansions. Compute interactions of each { y_k} in subinterval i of level nlevs with all {x_j} in subinterval i-1,i,i+1. Add all these interactions with the respective local expansions.  

				Complexity: O(3Ns)

			 
					
				
1				Total Complexity of FMM
				
				
    O(1+1+Np+(22Np^2/s)+ Np + (82Np^2/s)+Np+3Ns)
        =     O(2+3Np+3Ns+(10Np^2/s))
        =     O(2+3Np+6Np+(10Np^2/2p))
        =     O(2+9Np+(5Np))
        =     O(2+14Np)
        =     O(Np)
        =     O(N  log(1/ϵ))
    				, Where p = log(1/ϵ)

				
			
	
