WO2024009469A1 - 行列単純化装置、行列単純化方法、およびプログラム - Google Patents

行列単純化装置、行列単純化方法、およびプログラム Download PDF

Info

Publication number
WO2024009469A1
WO2024009469A1 PCT/JP2022/027010 JP2022027010W WO2024009469A1 WO 2024009469 A1 WO2024009469 A1 WO 2024009469A1 JP 2022027010 W JP2022027010 W JP 2022027010W WO 2024009469 A1 WO2024009469 A1 WO 2024009469A1
Authority
WO
WIPO (PCT)
Prior art keywords
matrix
norm
low
karyotype
weighted
Prior art date
Application number
PCT/JP2022/027010
Other languages
English (en)
French (fr)
Inventor
崇元 佐々木
幸浩 坂東
正樹 北原
Original Assignee
日本電信電話株式会社
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by 日本電信電話株式会社 filed Critical 日本電信電話株式会社
Priority to PCT/JP2022/027010 priority Critical patent/WO2024009469A1/ja
Publication of WO2024009469A1 publication Critical patent/WO2024009469A1/ja

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization

Definitions

  • the present invention relates to a matrix simplification device, a matrix simplification method, and a program.
  • low rank means that the rank (X) of the matrix X derived from data or physical phenomena is compared with the size of X (here, the smaller of the number of rows or the number of columns of matrix X). It has a small nature.
  • a mathematical programming problem is formulated to minimize the rank of the focused matrix, and a solution to this problem is found to achieve the desired analysis.
  • the function rank is discontinuous, non-differentiable, and non-convex, and the planning problem based on it is an NP-hard combinatorial optimization. For this reason, a relaxation approach that regularizes the karyotype norm is widely used instead of the rank function. Since the karyotype norm is a convex envelope of the rank function, minimizing the karyotype norm can indirectly improve low rank. Also, an approach that regularizes the weighted karyotype norm instead of the karyotype norm is also widely used. The weighted karyotype norm is more effective in increasing low rank than the karyotype norm.
  • karyotype norms and weighted karyotype norms are solved by iterative calculations of proximity separation methods such as Alternating Direction Method of Multipliers (ADMM), which are proximal mappings of karyotype norms and weighted karyotype norms.
  • ADMM Alternating Direction Method of Multipliers
  • SVT Singular value thresholding
  • the karyotype norm regularization problem formulated by each analysis method is classified into (1) a problem of regularizing a small number of large matrices, and (2) a problem of regularizing a large number of small matrices.
  • Applications of the former (1) include, for example, robust principal component analysis, missing value estimation/interpolation, optical flow estimation, dynamic MRI (Magnetic Resonance Imaging) analysis, and genome analysis.
  • the latter (2) is used for graph simplification, false color removal, and the like.
  • Patent Document 1 and Non-Patent Document 3 propose a fast parallel SVT (FPSVT) that calculates the SVT of the karyotype norm using data in parallel while suppressing the amount of calculation.
  • FPSVT fast parallel SVT
  • the karyotype norm is expressed by the L ⁇ ,2 mixed norm without using singular values, thereby realizing SVT calculation that does not require SVD and reducing the amount of calculation.
  • a parallelization algorithm that processes data in parallel can be easily implemented, and many matrices can be processed simultaneously.
  • Patent Document 1 Although the SVT of the karyotype norm can be determined at high speed, the SVT of the weighted karyotype norm cannot be determined. In order to further enhance low-rank properties in low-rank modeling, a technique is required to quickly obtain the SVT of the weighted karyotype norm.
  • the present invention was made based on the above-mentioned recognition of the problem, and allows SVT (singular value threshold processing) of weighted karyotype norms to be calculated in parallel with less calculation amount, and to reduce the matrix at higher speed.
  • the present invention aims to provide a matrix simplification device, a matrix simplification method, and a program that enable ranking.
  • a matrix simplification device includes a rank approximation unit.
  • the vectorization step and the weighted karyotype norm by the weight w of the matrix Y [w 1 , w 2 ] T , when the rotation matrix on the image ImY is R, (y 1 - Ry 2 , y 1 + Ry 2 ) Low-rank approximation that performs low-rank approximation of the matrix Y by calculating the SVT (Singular Value Thresholding) of the weighted karyotype norm based on the group cusp quasi-norm based on
  • the rotation matrix on the image ImY is R, (y 1 - Ry 2 , y 1 + Ry 2 ) and w, and by calculating the SVT (Singular Value Thresholding) of the weighted karyotype norm, a low-rank approximation of the matrix Y is performed.
  • This is a program for functioning as a rank approximation unit.
  • the SVT single value threshold processing
  • FIG. 2 is a diagram illustrating an embodiment of the present invention, and is a graph showing the relationship between vectors y 1 , y 2 , Ry 2 , and ⁇ Ry 2 in space ImY.
  • FIG. 2 is a schematic diagram showing Algorithm 1 (SVT calculation method when M ⁇ 3).
  • FIG. 2 is a block diagram showing an example of a schematic functional configuration of a matrix simplification device.
  • FIG. 2 is a schematic diagram showing an example of the structure of main data processed by the matrix simplification device.
  • FIG. 2 is a schematic diagram showing an example of the structure of main data processed by the matrix simplification device.
  • FIG. 2 is a schematic diagram showing an example of the structure of main data processed by the matrix simplification device.
  • FIG. 2 is a schematic diagram showing an example of the structure of main data processed by the matrix simplification device.
  • 3 is a flowchart illustrating an example of a
  • the weighted karyotype norm is defined as the sum of the weights w i of the singular values ⁇ i (X). That is, as shown in the following equation (1).
  • the present embodiment deals with an optimization problem in which a large number of small matrices are regularized with a weighted karyotype norm. This problem has the form of equation (2) below.
  • the second term of the objective function is a low-rank regularization function based on the sum of weighted karyotype norms, and the first term is another regularization or fidelity function.
  • the linear mapping L i is a linear mapping for generating a matrix that regulates low rank.
  • Non-Patent Document 4 and Non-Patent Document 5 use a proximate separation method such as ADMM or Primal Dual Splitting because the function f is differentiable or proximate mapping calculation is possible.
  • Non-patent Document 4 Takamoto Sasaki, Ryuichi Tanida, Jun Shimizu, “Graph shape simplification by local linear approximation of graph signals,” 15th Information Science and Technology Forum (FIT2016), Volume 3, pp.1-4 ,Sept.2016.
  • Non-patent Document 5 S. Ono and I. Yamada, “Color-line regularization for color artifact removal,” IEEE Transactions on Computational Imaging, vol.2, no.3, pp.204-217, 2016.
  • proximity mapping computation possible is as follows. That is, the proximity mapping prox g at the point (vector) y for the function g is defined as in the following equation (3). Efficient methods of calculating proximity maps are known for some functions including weighted karyotype norms, and such functions are called proximity mapping computable.
  • the above proximity separation method is an iterative calculation procedure that generates a convergence sequence to the optimal solution, and in each iteration, the solution is updated based on primary information such as the gradient of the function and the proximity mapping, and the objective function is minimized.
  • the proximity mapping of g(Y 1 , . . . , Y L ) can be calculated by independently calculating the proximity mapping of each element constituting the right side, and the following equation (5) holds true.
  • Each element on the right side is a proximity map of the weighted karyotype norm and is equal to the SVT of the weighted karyotype norm.
  • the weighted karyotype norm is expressed as
  • the SVT of the weighted karyotype norm is calculated by the following equation (6), where Y is the input matrix and w is the weight vector.
  • the function SVD( ⁇ ) is singular value decomposition (SVD).
  • (.) + is a ramp function that clips each element of the input to a non-negative value.
  • diag( ⁇ ) is a function that generates a diagonal matrix from an input vector.
  • the proximity mapping of the function g(Y 1 , Y 2 , . . . , Y L ) can be calculated by L times of SVT, and each SVT requires SVD calculation using a small matrix as input.
  • iterative calculations are performed using the proximity separation method, most of the calculation time is spent on SVD calculation, which requires a large amount of calculation.
  • each of the matrices Y 1 , Y 2 , . . . , Y L can be processed in L parallel. Note that when the size of each matrix is small in the problem being handled, the overhead of processing allocation and memory loading is relatively large, and the improvement effect of parallelization may be limited.
  • FPSVT Fast Parallel SVT
  • This method is derived from the feature that the weighted karyotype norm is expressed by the group cusp quasi-norm. This feature allows the weighted karyotype norm to be expressed without using singular values.
  • SVT can be expressed without using SVD. Since most of this SVT can be described by linear transformation, data-parallel algorithms can be derived.
  • the SVT calculation method obtained in this embodiment reduces the amount of calculation and can be executed as a data parallel algorithm.
  • the singular value sum ⁇ 1 + ⁇ 2 and the singular value difference ⁇ 1 ⁇ 2 are as expressed by the following equation (8) according to Non-Patent Document 3.
  • R is a rotation matrix, and this rotation matrix R rotates the vector on the image ImY by ⁇ /2 [rad] around the origin along ImY.
  • the positive and negative directions of rotation as shown in Figure 1, the direction that takes you from y 1 to y 2 in the shortest distance is the positive direction (first direction), and the opposite direction is the negative direction (second direction). ).
  • a and ⁇ are expressed by the following equations (10-1) and (10-2).
  • Non-Patent Document 6 Takamoto Sasaki, Yukihiro Bando, Masaki Kitahara, “Edge-preserving smoothing that allows simultaneous control of gradation and edges based on cusp star total variation regularization,” Proceedings of the 36th Signal Processing Symposium, Electronic Information Telecommunications Society Signal Processing Research Committee, 2021, pp.181-186.
  • prox ⁇ (y 1 ', y 2 ') is It is.
  • sgn( ⁇ ) is a sign function
  • ⁇ points are element products
  • ⁇ is a vector obtained by sorting [
  • proj ⁇ m ( ⁇ ) is a projection onto a non-increasing monotone cone and can be solved using the Pool Adjacent Violators Algorithm with a computational complexity of O(N).
  • the matrix R can be calculated as shown in the following equation (22).
  • the SVT of the weighted karyotype norm can be determined by the procedure of determining
  • (i), (ii), and (v) can be calculated using product-sum calculations, (iii) can be calculated using four arithmetic calculations and radical calculations, and (iv) can be calculated using four arithmetic calculations and clipping. Therefore, it is an SVT of many input matrices Y 1 ,...,Y L can be obtained in parallel at high speed by SIMD calculation.
  • FIG. 2 is a schematic diagram showing an algorithm for calculating SVT (Singular Value Thresholding) when M ⁇ 3. This algorithm is written in pseudo code. The algorithm will be explained below with reference to this diagram.
  • the inputs are an M ⁇ 2 (M rows and 2 columns) matrix Y and an SVT weight vector w.
  • the column vectors (M dimension) of the first and second columns of Y are expressed as y 1 and y 2 , respectively.
  • w [w 1 , w 2 ] T satisfies w 2 ⁇ w 1 ⁇ 0.
  • the output is an SVT based on the weight w of the matrix Y.
  • the output matrix is denoted by Z.
  • assignments to variables d, e, and f are made.
  • the value ac-b 2 is assigned to the variable d
  • the value ⁇ d is assigned to the variable e
  • the value a+c is assigned to the variable f.
  • the value of the variable d is the determinant of the matrix Y T Y
  • the value of the variable e is its square root
  • the value of the variable f is the trace of the matrix Y T Y.
  • assignments to variables ⁇ 1 and ⁇ 2 are performed.
  • the value (g+h)/2 is assigned to the variable ⁇ 1
  • the value (gh ⁇ h)/2 is assigned to the variable ⁇ 2 .
  • the value of the variable ⁇ 1 is the maximum singular value of the matrix Y
  • the value of the variable ⁇ 2 is the minimum singular value of the matrix Y.
  • the inputs are a 2 ⁇ 2 (2 rows and 2 columns) matrix Y and the weight W of the SVT.
  • the column vectors (two-dimensional) of the first and second columns of Y are expressed as y 1 and y 2, respectively, and each element of the matrix Y is given a suffix (subscript) in which the row number and column number are arranged in this order.
  • w [w 1 , w 2 ] T satisfies w 2 ⁇ w 1 ⁇ 0.
  • the output is an SVT based on the weight w of the matrix Y.
  • the output matrix is denoted by Z.
  • the value y 1,1 y 2,2 -y 1,2 y 2,1 is assigned to the variable d
  • is assigned to the variable e
  • the value a+c is assigned to the variable f.
  • the value of the variable d is the determinant of the matrix Y
  • the value of the variable e is its absolute value and the square root of the determinant of the matrix YTY
  • the value of the variable f is the trace of the matrix YTY .
  • assignments to variables ⁇ 1 and ⁇ 2 are performed.
  • the value (g+h)/2 is assigned to the variable ⁇ 1
  • the value (gh ⁇ h)/2 is assigned to the variable ⁇ 2 .
  • the value of the variable ⁇ 1 is the maximum singular value of the matrix Y
  • the value of the variable ⁇ 2 is the minimum singular value of the matrix Y.
  • the amount of calculation for SVT calculation using the algorithm of this embodiment will be considered.
  • the SVT calculation method (prior art) using SVD requires 24M+160 floating point operations to obtain the SVD, 2 times for threshold processing, and 6M+4 floating point operations to obtain the matrix product. That is, a total of 30M+166 floating point operations are required.
  • Algorithm 1 and Algorithm 2 can be efficiently processed in parallel using a parallel architecture such as SIMD.
  • FIG. 4 is a block diagram showing a schematic functional configuration of the matrix simplification device 1 according to this embodiment.
  • the matrix simplification device 1 includes an input section 11, a vectorization section 12, a low-rank approximation section 13, and an output section 14.
  • each unit includes a storage unit therein for storing data as necessary.
  • This storage unit is realized using storage means such as a semiconductor memory or a magnetic hard disk device.
  • the functions of each part may be realized by a computer and a program.
  • the input unit 11 acquires data in a matrix of M rows and 2 columns or 2 rows and M columns from the outside.
  • M is an integer of 2 or more.
  • the elements of this matrix are numerical values (scalar values).
  • the input unit 11 transposes the input matrix as necessary. In other words, when the vectorization unit 12 and the low-rank approximation unit 13 in the subsequent stage are configured to process only a matrix of either M rows and 2 columns or 2 rows and M columns, and the input When the input matrix does not match the format (that is, the length and width are reversed), the input unit 11 transposes the input matrix. This allows the matrix simplification device 1 to process either a matrix with M rows and 2 columns or 2 rows and M columns.
  • the vectorization unit 12 and the low-rank approximation unit 13 will be described as processing a matrix with M rows and 2 columns. However, this may be for processing a matrix of 2 rows and M columns, and the essential processing content remains the same.
  • the vectorization unit 12 vectorizes a matrix (denoted as Y) with M rows and 2 columns.
  • the vectorization unit 12 passes these vectors y 1 and y 2 to the low-rank approximation unit 13 .
  • the vectorization unit 12 extracts an M-dimensional first vector and a second vector corresponding to each column from an input matrix of M rows and 2 columns (or extracts an M-dimensional first vector and a second vector corresponding to each row from an input matrix of 2 rows and M columns). ).
  • FPSVT Fast Parallel SVT
  • the output unit 14 outputs the data of the matrix Z obtained by the low-rank approximation unit 13 to the outside. That is, the matrix simplification device 1 performs low-rank simplification on the input matrix Y, and outputs the matrix Z that is the result. Note that when the input unit 11 transposes the matrix, the output unit 14 transposes the matrix Z obtained by the low-rank approximation unit 13 again and then outputs it. This allows the size of the input matrix (number of rows and columns) to match the size of the output matrix.
  • FIG. 5A to 5C are schematic diagrams showing the structure of main data processed by the matrix simplification device 1.
  • FIG. 5A shows the data structure of the input matrix Y.
  • a matrix with M rows and 2 columns is shown.
  • FIG. 5B shows the data structure of vectors y 1 and y 2 vectorized based on matrix Y.
  • Vectors y 1 and y 2 are each M-dimensional column vectors.
  • the first column of matrix Y corresponds to vector y1
  • the second column corresponds to vector y2 .
  • FIG. 5C shows the data structure of the output matrix Z.
  • the matrix Z also shows a matrix with M rows and 2 columns.
  • FIG. 6 is a flowchart showing the process flow of the matrix simplification device 1 according to this embodiment.
  • the input unit 11 transposes the input matrix to have M rows and 2 columns (step S102).
  • the vectorization unit 12 divides the M-row, 2-column matrix Y into two M-dimensional column vectors y 1 (first vector) and y 2 (second vector) (step S103).
  • the low-rank approximation unit 13 determines whether M ⁇ 3 (step S104). If Yes (M ⁇ 3) in step S104, the low-rank approximation unit 13 calculates y 1 T y 1 , y 1 T y 2 , y 2 T y 2 and sets the variables a, b, c, respectively. (Step S105). These a, b, and c are all scalar values.
  • the low-rank approximation unit 13 calculates ac-b 2 , ⁇ d, and a+c, and assigns them to variables d, e, and f, respectively (step S106).
  • the value of the variable d is the determinant of the matrix Y T Y
  • the value of the variable e is its square root
  • the value of the variable f is the trace of the matrix Y T Y.
  • the low-rank approximation unit 13 calculates ⁇ (f+2e) and ⁇ (f ⁇ 2e) and assigns them to variables g and h, respectively (step S107).
  • the value of the variable g is the sum of the singular values of the matrix Y
  • the value of the variable h is the difference of the singular values of the matrix Y.
  • the low-rank approximation unit 13 calculates (g+h)/2 and (gh)/2 and assigns them to variables ⁇ 1 and ⁇ 2 , respectively (step S108).
  • the value of the variable ⁇ 1 is the maximum singular value of the matrix Y
  • the value of the variable ⁇ 2 is the minimum singular value of the matrix Y.
  • the low-rank approximation unit 13 calculates that if g is not 0, is calculated and assigned to the variable ⁇ . On the other hand, if g is 0, the low-rank approximation unit 13 assigns 0 to the variable ⁇ (step S109).
  • the low-rank approximation unit 13 calculates that if h is not 0, is calculated and assigned to the variable ⁇ . On the other hand, if h is 0, the low-rank approximation unit 13 assigns 0 to the variable ⁇ (step S110).
  • the low rank approximation unit 13 is calculated and assigned to variable Z (step S111).
  • step S104 the low-rank approximation unit 13 calculates y 1 T y 1 and y 2 T y 2 and assigns them to variables a and c, respectively (step S112 ). Both a and c are scalar values.
  • the low-rank approximation unit 13 calculates y 1,1 y 2,2 ⁇ y 1,2 y 2,1 ,
  • the value of the variable d is the determinant of the matrix Y
  • the value of the variable e is its absolute value and the square root of the determinant of the matrix YTY
  • the value of the variable f is the trace of the matrix YTY .
  • the low-rank approximation unit 13 calculates ⁇ (f+2e) and ⁇ (f ⁇ 2e) and assigns them to variables g and h, respectively (step S114).
  • the value of the variable g is the sum of the singular values of the matrix Y
  • the value of the variable h is the difference of the singular values of the matrix Y.
  • the low-rank approximation unit 13 calculates (g+h)/2 and (gh)/2 and assigns them to variables ⁇ 1 and ⁇ 2 , respectively (step S115).
  • the value of the variable ⁇ 1 is the maximum singular value of the matrix Y
  • the value of the variable ⁇ 2 is the minimum singular value of the matrix Y.
  • the low-rank approximation unit 13 calculates that if g is not 0, is calculated and assigned to the variable ⁇ . On the other hand, if g is 0, the low-rank approximation unit 13 assigns 0 to the variable ⁇ (step S116).
  • the low-rank approximation unit 13 calculates that if h is not 0, is calculated and assigned to the variable ⁇ . On the other hand, if h is 0, the low-rank approximation unit 13 assigns 0 to the variable ⁇ (step S117).
  • the low rank approximation unit 13 is calculated and assigned to variable Z (step S118).
  • step S111 or step S118 the output unit 14 outputs the value of the variable Z (Z T if the input matrix has been transposed in step S102). This concludes the explanation of the process flow of the matrix simplification device 1 in FIG. 6.
  • All or part of the functions of the matrix simplification device 1 in each of the embodiments described above may be realized by a computer.
  • a program for realizing this function may be recorded on a computer-readable recording medium, and the program recorded on the recording medium may be read into a computer system and executed.
  • the "computer system” herein includes hardware such as an OS and peripheral devices.
  • “computer-readable recording media” refers to portable media such as flexible disks, magneto-optical disks, ROM, CD-ROM, DVD-ROM, and USB memory, and storage devices such as hard disks built into computer systems. Say something.
  • a "computer-readable recording medium” refers to a storage medium that dynamically stores a program for a short period of time, such as a communication line when transmitting a program via a network such as the Internet or a communication line such as a telephone line. In that case, it may also include something that retains a program for a certain period of time, such as a volatile memory inside a computer system that is a server or client. Further, the program may be one for realizing a part of the above-mentioned functions, or may be one that can realize the above-mentioned functions in combination with a program already recorded in the computer system.
  • the matrix simplification device 1 is capable of executing a high-speed calculation algorithm of singular value threshold processing (SVT). This makes it possible to quickly find a solution to an optimization problem that involves low-rank regularization of many small matrices.
  • the matrix simplification device 1 can perform rank reduction processing at high speed with a small amount of calculation for a matrix with M rows and 2 columns or 2 rows and M columns (M ⁇ 2). Become.
  • the matrix simplification device 1 can execute rank reduction processing in parallel for a plurality of matrices. That is, the matrix simplification device 1 realizes an FPSVT (Fast Parallel SVT) algorithm for a matrix with M rows and 2 columns or 2 rows and M columns.
  • FPSVT Fast Parallel SVT
  • the matrix simplifier 1 has realized a data-parallel SVT calculation method that does not require singular value decomposition (SVD). Furthermore, as a result of evaluation experiments using actual data, it was confirmed that SVT could be calculated up to 77 times faster while having higher calculation accuracy than conventional methods.
  • SVD singular value decomposition
  • the present invention can be applied to techniques that handle optimization problems in which a large number of small matrices are regularized with a weighted karyotype norm. Furthermore, the present invention can be applied to techniques for performing data analysis based on low-rank characteristics latent in data and physical phenomena, such as graph simplification and false color removal.

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Computational Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • Computing Systems (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Image Processing (AREA)

Abstract

行列単純化装置は、Ｍ行２列行列Ｙ＝［ｙ１，ｙ２］又は２行Ｍ列行列ＹＴ＝［ｙ１，ｙ２］Ｔを、２つのＭ次元ベクトルｙ１，ｙ２とするベクトル化部と、重みｗ＝［ｗ１，ｗ２］Ｔによる行列Ｙの加重核型ノルムが、像ＩｍＹ上の回転行列をＲとするとき、［ｙ１－Ｒｙ２，ｙ１＋Ｒｙ２］と前記重みｗに基づくグループ尖星準ノルムで表されることに基づいて、前記加重核型ノルムのＳＶＴ（Singular Value Thresholding）を算出することにより、行列Ｙの低ランク近似化を行う低ランク近似化部と、を備える。

Description

行列単純化装置、行列単純化方法、およびプログラム
　本発明は、行列単純化装置、行列単純化方法、およびプログラムに関する。
　データや物理現象に潜在する低ランク性に基づくデータ解析手法（以降、低ランクモデリング）が、コンピュータビジョン、画像処理、ゲノムデータ解析などの多くの分野で近年活発に研究されている。ここで、低ランク性とは、データや物理現象から導かれる行列Ｘについて、その階数ｒａｎｋ（Ｘ）がＸのサイズ（ここでは、行列Ｘの行数または列数のうちの小さい方）と比較して小さい性質である。低ランクモデリングでは、着目した行列の階数（rank）を最小化する数理計画問題を立式し、この問題の解を求めて所望の解析を実現する。
　関数rankは不連続で微分不能、非凸であり、これに基づく計画問題はＮＰ困難な組合せ最適化となる。このため、rank関数の代わりに核型ノルムを正則化する緩和アプローチが広く用いられる。核型ノルムはrank関数の凸包絡であるため、核型ノルムを最小化することで間接的に低ランク性を高められる。また核型ノルムの代わりに加重核型ノルムを正則化するアプローチも広く用いられる。加重核型ノルムは低ランク性を高める効果が核型ノルムより高い。
　核型ノルムや加重核型ノルムを含む最適化問題は、Alternating Direction Method of Multipliers（ＡＤＭＭ）等の近接分離法の反復計算で解を求められ、核型ノルムや加重核型ノルムの近接写像である特異値閾値処理（Singular Value Thresholding；ＳＶＴ）が繰り返し実行される。
　しかし、ＳＶＴ算出には計算量の大きい特異値分解（Singular Value Decomposition；ＳＶＤ）の算出が必要なため、解析結果を得るのに多くの計算時間を要する。ここで、各解析法で立式される核型ノルム正則化問題を、（１）少数の大型行列を正則化する問題と、（２）多数の小型行列を正則化する問題に分類する。前者（１）の用途は、例えば、ロバスト主成分分析や、欠損値推定・補間や、オプティカルフロー推定や、ダイナミックＭＲＩ（Magnetic Resonance Imaging）解析や、ゲノム解析等である。また、後者（２）の用途は、グラフ単純化や、偽色除去等である。
　前者（１）の問題の場合には、計算量を抑えてＳＶＴを高速化する手法がいくつか提案されている。J.F.Caiらは、行列を事前に完全直交分解（Complete Orthogonal Decomposition；ＣＯＤ）した後にニュートン法で反復更新し、ＳＶＤを行わずにＳＶＴを求める高速ＳＶＴ（Fast SVT；ＦＳＶＴ）を提案している（非特許文献１）。またT.H.Ohらは、大型行列を直交行列と小型のコア行列の積に近似することで、ＳＶＤの入力サイズを小さくして高速化する高速ランダム化ＳＶＴ（Fast Randomized SVT；ＦＲＳＶＴ）を提案している（非特許文献２）。ここに挙げるいずれの手法も、入力行列のサイズが大きい（行数および列数がそれぞれ５００～２０００程度）ときに計算量を抑え、大幅な速度改善効果を示す。
　一方で後者（２）の問題の場合、上記手法による高速化の効果は限定的であると考えられる。ＦＳＶＴを用いると、入力が小型の場合、直接のＳＶＤ計算と比較してＣＯＤとニュートン法の計算量が大きいという問題がある。また、ＦＲＳＶＴを用いると、入力が小型の場合、コア行列の縮退効果が小さいために高速化できず、また近似法であるため計算誤差が大きいという問題がある。加えて、これらの手法は多数の行列を同時に処理するデータ並列のアプローチを取れず、近年の並列アーキテクチャの計算資源を有効に活用できないという問題もある。
　特許文献１および非特許文献３では、計算量を抑えながらデータを並列に核型ノルムのＳＶＴを算出する高速並列ＳＶＴ（Fast Parallel SVT；ＦＰＳＶＴ）が提案されている。ＦＰＳＶＴでは、特異値を用いずにＬ∞，２混合ノルムで核型ノルムを表現することで、ＳＶＤが不要なＳＶＴ計算を実現し、計算量を削減している。加えて、データを並列に処理する並列化アルゴリズムを容易に実現でき、多数の行列について同時処理が可能である。
特許第６８１０００３号公報
J.-F. Cai, O. Stanley, "Fast singular value thresholding without singular value decomposition," Methods and Applications of Analysis, vol.20, no.4, pp.335-352, Dec. 2013. T. H. Oh, Y. Matsushita, Y. W. Tai, and I. S. Kweon, "Fast randomized singular value thresholding for nuclear norm minimization," 2015 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pp.4484-4493, June 2015. 佐々木崇元, 北原正樹, 清水淳, "低ランク最適化のための高速特異値閾値処理の数理,"第16回 情報科学技術フォーラム (FIT), 第1分冊, 2017, pp.5－12.
　しかしながら特許文献１、非特許文献３の技術では、核型ノルムのＳＶＴを高速に求められるが、加重核型ノルムのＳＶＴを求められない。低ランクモデリングにおいて低ランク性をより高めるために、加重核型ノルムのＳＶＴを高速に求める技術が求められている。
　本発明は、上記の課題認識に基づいて行われたものであり、加重核型ノルムのＳＶＴ（特異値閾値処理）についてより少ない計算量でかつ並列に求めることができ、より高速に行列の低ランク化をすることが可能となる行列単純化装置、行列単純化方法、およびプログラムを提供することを目的とする。
　本発明の一態様は、Ｍ行２列行列Ｙ＝［ｙ１，ｙ２］又は２行Ｍ列行列ＹＴ＝［ｙ１，ｙ２］Ｔを、２つのＭ次元ベクトルｙ１，ｙ２とするベクトル化部と、重みｗ＝［ｗ１，ｗ２］Ｔによる行列Ｙの加重核型ノルムが、像ＩｍＹ上の回転行列をＲとするとき、［ｙ１－Ｒｙ２，ｙ１＋Ｒｙ２］と前記重みｗに基づくグループ尖星準ノルムで表されることに基づいて、前記加重核型ノルムのＳＶＴ（Singular Value Thresholding）を算出することにより、行列Ｙの低ランク近似化を行う低ランク近似化部と、を備える行列単純化装置である。
　本発明の一態様は、Ｍ行２列行列Ｙ＝［ｙ１，ｙ２］又は２行Ｍ列行列ＹＴ＝［ｙ１，ｙ２］Ｔを、２つのＭ次元ベクトルｙ１，ｙ２とするベクトル化ステップと、行列Ｙの重みｗ＝［ｗ１，ｗ２］Ｔによる加重核型ノルムが、像ＩｍＹ上の回転行列をＲとするとき、（ｙ１－Ｒｙ２，ｙ１＋Ｒｙ２）とｗに基づくグループ尖星準ノルムで表されることに基づいて、前記加重核型ノルムのＳＶＴ（Singular Value Thresholding）を算出することにより、行列Ｙの低ランク近似化を行う低ランク近似化ステップと、を含む行列単純化方法である。
　本発明の一態様は、コンピュータを、Ｍ行２列行列Ｙ＝［ｙ１，ｙ２］又は２行Ｍ列行列ＹＴ＝［ｙ１，ｙ２］Ｔを、２つのＭ次元ベクトルｙ１，ｙ２とするベクトル化部、行列Ｙの重みｗ＝［ｗ１，ｗ２］Ｔによる加重核型ノルムが、像ＩｍＹ上の回転行列をＲとするとき、（ｙ１－Ｒｙ２，ｙ１＋Ｒｙ２）とｗに基づくグループ尖星準ノルムで表されることに基づいて、前記加重核型ノルムのＳＶＴ（Singular Value Thresholding）を算出することにより、行列Ｙの低ランク近似化を行う低ランク近似化部、として機能させるためのプログラムである。
　本発明によれば、加重核型ノルムのＳＶＴ（特異値閾値処理）についてより少ない計算量でかつ並列に求めることができ、より高速に行列の低ランク化をすることが可能となる。
本発明の実施形態を説明する図であり、空間ＩｍＹにおけるベクトルｙ１，ｙ２，Ｒｙ２，－Ｒｙ２の関係を示すグラフである。 アルゴリズム１（Ｍ≧３の場合のＳＶＴ算出法）を示す概略図である。 アルゴリズム２（Ｍ＝２の場合のＳＶＴ算出法）を示す概略図である。 行列単純化装置の概略機能構成の例を示すブロック図である。 行列単純化装置が処理する主要なデータの構成の例を示す概略図である。 行列単純化装置が処理する主要なデータの構成の例を示す概略図である。 行列単純化装置が処理する主要なデータの構成の例を示す概略図である。 行列単純化装置の処理の流れの例を示すフローチャートである。
　次に、本発明の実施形態について、図面を参照しながら説明する。以下において、実施形態に共通する低ランクモデリングについて説明し、その後で行列単純化装置１等の具体的な構成を説明する。
［１．加重核型ノルム正則化に基づく低ランクモデリング］
　本実施形態が対象とする加重核型ノルム正則化問題は、次に説明する通りの問題である。
　行列Ｘ∈ＲＭ×Ｎ（ここで、Ｒは実数の集合。ＲＭ×Ｎは実数を要素とするＭ行Ｎ列の行列の集合。以下においても同様。）の特異値をσｉ（Ｘ）、ｉ＝１，２，・・・，Ｋとする。ただし、Ｋ＝ｍｉｎ（Ｍ，Ｎ）であり、ｍｉｎ（）は引数の最小値を返す関数である。なお、以後において、Ｍ行Ｎ列の行列のサイズを単に「Ｍ×Ｎ」と表す場合がある。
　このとき、加重核型ノルムは特異値σｉ（Ｘ）の重み付けｗｉの和として定義される。即ち、次の式(1)の通りである。
Figure JPOXMLDOC01-appb-M000008
 この加重核型ノルムを正則化することで低ランク性を推進できる。
　本実施形態は、多数の小型行列が加重核型ノルムで正則化される最適化問題を取り扱う。この問題は次の式(2)の形式を有する。
Figure JPOXMLDOC01-appb-M000009
　式(2)で目的関数の第２項は加重核型ノルムの和による低ランク性の正則化関数であり、第１項はその他の正則化や忠実化関数である。線形写像Ｌｉは低ランク性を規約する行列を生成するための線形写像である。
　背景技術において述べた「（２）多数の小型行列を正則化する問題」を解くためには、関数ｆと線形写像Ｌｉに応じたアルゴリズムを選択する。そのアルゴリズムとして、非特許文献４や、非特許文献５では、関数ｆが微分可能あるいは近接写像計算可能であるため、ＡＤＭＭやPrimal Dual Splitting等の近接分離法を用いている。
［非特許文献４］佐々木崇元，谷田隆一，清水淳，“グラフ信号の局所線形近似によるグラフ形状単純化，” 第１５回情報科学技術フォーラム(FIT2016)，第３分冊，pp.1-4，Sept.2016.
［非特許文献５］S. Ono and I. Yamada, “Color-line regularization for color artifact removal,” IEEE Transactions on Computational Imaging, vol.2, no.3, pp.204-217, 2016.
　なお、ここで「近接写像計算可能」の意味は次の通りである。即ち、関数ｇに対し、点（ベクトル）ｙでの近接写像ｐｒｏｘｇは、次の式(3)の様に定義される。
Figure JPOXMLDOC01-appb-M000010
 　加重核型ノルムを含むいくつかの関数については近接写像の効率的な計算法が知られており、そのような関数を近接写像計算可能であると呼ぶ。
　上記の近接分離法は、最適解への収束列を生成する反復計算手続きで、各反復では関数の勾配や近接写像などの１次情報に基づき解を更新し、目的関数を最小化する。
　「多数の小型行列を正則化する問題」の加重核型ノルムの和の正則化関数については、補助変数Ｙｉ＝ＬｉＸ（ｉ＝１，…，Ｌ）を導入し、次の式(4)で表されるｇ（Ｙ１，…，ＹＬ）の近接写像を計算すれば良い。
Figure JPOXMLDOC01-appb-M000011
　ｇ（Ｙ１，…，ＹＬ）の近接写像は、右辺を構成する各要素の近接写像を独立に計算すれば良く、次の式(5)が成立する。
Figure JPOXMLDOC01-appb-M000012
 　右辺の各要素は加重核型ノルムの近接写像で、加重核型ノルムのＳＶＴに等しい。加重核型ノルムは上の式において｜｜・｜｜＊，ｗと表されている。
　加重核型ノルムのＳＶＴは、入力行列をＹ、重みベクトルをｗとして、次の式(6)で算出される。
Figure JPOXMLDOC01-appb-M000013
 　ここで、関数ＳＶＤ（・）は、特異値分解（ＳＶＤ）である。また、（・）＋は、入力の各要素を非負値にクリッピングするランプ関数である。ｄｉａｇ（・）は入力ベクトルから対角行列を生成する関数である。
　なお、本実施形態での特異値分解は、「thin SVD」である。即ち、Ｙ＝ＵΣＶＴと分解した時（Ｙ∈ＲＭ×Ｎ）、行列Ｕ，ＶはＫ＝ｍｉｎ（Ｍ，Ｎ）個の正規直交なベクトルから成り、（Ｋ＋１）番目以降の特異ベクトルの算出は省略する（Ｕ∈ＲＭ×Ｋ，Ｖ∈ＲＮ×Ｋ）。また、行列ΣはＫ×Ｋ対角行列である。数値解析ソフトＭＡＴＬＡＢ（MathWorks社，米国）で記述する場合は、「[U,S,V]=svd(Y,'econ');」である。
　即ち、関数ｇ（Ｙ１，Ｙ２，・・・，ＹＬ）の近接写像は、Ｌ回のＳＶＴにより算出でき、それぞれのＳＶＴでは小型行列を入力とするＳＶＤ算出が必要である。近接分離法で反復演算する際、計算時間の多くは計算量の大きいＳＶＤ算出に費やされる。
　さて、関数ｇ（Ｙ１，Ｙ２，・・・，ＹＬ）の近接写像はＬ個の行列Ｙ１，Ｙ２，・・・，ＹＬの各々の間における依存関係が無いため、タスク並列処理が可能である。つまり、行列Ｙ１，Ｙ２，・・・，ＹＬの各々についてＬ並列での処理が可能である。なお、扱う問題において各行列のサイズが小さい場合には、処理割当やメモリロードのオーバーヘッドが相対的に大きく、並列化による改善の効果が限定的である場合もある。
　そこで、オーバーヘッドが少ない、データ並列処理について説明する。アルゴリズムをデータ並列化できれば、Single Instruction Multiple Data（ＳＩＭＤ）等のデータ並列アーキテクチャを用いた実装により処理を高速化できる。しかし、ＳＶＤの算出は逐次的であり、且つ行列要素の参照位置や処理内容が入力Ｙｉに依存するため、行列間で共通に計算できる処理が少なく、データ並列化は本質的に困難である。
　以上より、多数の小型行列を低ランクに正則化するモデルは最適解を高速に得るのが困難で、その原因は、ＳＶＤの多大な計算量と並列化の難しさにある。
［２．Fast Parallel SVT（ＦＰＳＶＴ）］
　ここでは、多数のＳＶＴ計算を高速に算出するＦＰＳＶＴを導出する。この手法は加重核型ノルムがグループ尖星準ノルムにより表現されるという特徴に基づいて導かれる。この特徴により、特異値を用いずに加重核型ノルムを表現できる。またＳＶＤを用いずにＳＶＴを表現できる。そしてこのＳＶＴは、ほとんどが線形変換で記述できるため、データ並列なアルゴリズムを導ける。本実施形態で得られるＳＶＴ算出法では計算量が削減されており、かつデータ並列なアルゴリズムとして実行可能である。
　入力行列のサイズをＭ×２とする次の式(7)が成り立つことから、入力が２×Ｎのサイズの行列のＳＶＴの算出は、Ｎ×２のサイズの行列のＳＶＴ算出の前後に、式(7)の転置処理を施すことで実現できる。
Figure JPOXMLDOC01-appb-M000014
　まず特異値和σ１＋σ２と特異値差σ１－σ２は、非特許文献３により次の式(8)で表す通りである。なお式(8)において複号同順である。
Figure JPOXMLDOC01-appb-M000015
 ここでＲは回転行列であり、この回転行列Ｒは、像ＩｍＹ上のベクトルを、原点回りにＩｍＹに沿ってπ／２［ｒａｄ］回転させる。ただし回転方向の正負については、図１に示すように、ｙ１からｙ２に最短で辿り着く方向を回転の正方向（第１の方向）とし、その反対方向を負方向（第２の方向）とする。
　なお、上記の像ＩｍＹについて、次の通りである。即ち、行列Ａ∈ＲＭ×Ｎに対し、部分空間ＩｍＡ＝｛Ａｘ｜ｘ∈ＲＮ｝⊂ＲＭをＡの像という。
　これより、加重核型ノルムについて次の式(9)が成り立つ。
Figure JPOXMLDOC01-appb-M000016
ここで、Ａとωは次の式(10-1)，(10-2)で表される。
Figure JPOXMLDOC01-appb-M000017
 また式(9)において関数vecは入力行列Ｙ＝［ｙ１，ｙ２］∈ＲＭ×２を並べ替えて列ベクトル［ｙ１ Ｔ，ｙ２ Ｔ］Ｔ∈Ｒ２Ｍを出力する線形変換である。
　グループ尖星準ノルムΓω（ｙ’１，ｙ’２）は（非特許文献６）、次の式(11)で表される。
Figure JPOXMLDOC01-appb-M000018
 ここで｜｜ｙ’｜｜［１］と｜｜ｙ’｜｜［２］は、次の式(12)で表される。
Figure JPOXMLDOC01-appb-M000019
［非特許文献６］佐々木崇元, 坂東幸浩, 北原正樹, “尖星全変動正則化に基づくグラデーションとエッジ同時制御可能なエッジ保存平滑化,” 第36回信号処理シンポジウム講演論文集, 電子情報通信学会信号処理研究専門委員会, 2021, pp.181－186.
　これより近接写像は、
Figure JPOXMLDOC01-appb-M000020
 である。
　またｐｒｏｘΓω（ｙ１’，ｙ２’）は、
Figure JPOXMLDOC01-appb-M000021
 である。
ここで
Figure JPOXMLDOC01-appb-M000022
 である。ｓｇｎ（・）は符号関数、○に点は要素積、｜ｘ｜↓は［｜ｘ１｜，…，｜ｘN｜］Ｔを降順にソートしたベクトルである。ｐｒｏｊκｍ（・）は非増加単調錘への射影でPool Adjacent Violators AlgorithmによってＯ（Ｎ）の計算量で解ける。
ｐｒｏｘΩω（ａ，ｂ）において、ａ≧ｂ≧０，ω１＋ω２≧０，ω２≦０を考えると、ｓｇｎ関数、ソート、絶対値は必要がなくなり、
Figure JPOXMLDOC01-appb-M000023
 である。
これを計算すると、
Figure JPOXMLDOC01-appb-M000024
 である。
ゆえに
Figure JPOXMLDOC01-appb-M000025
 である。
よって、
Figure JPOXMLDOC01-appb-M000026
 である。
ここでα，βは
Figure JPOXMLDOC01-appb-M000027
 である。
　従って、
Figure JPOXMLDOC01-appb-M000028
 が成立する。
　行列Ｒは非特許文献３によれば、次の式(22)の通り計算できる。
Figure JPOXMLDOC01-appb-M000029
　従って、Ｍ≧３のときのＲ［ｙ２，－ｙ１］は、
Figure JPOXMLDOC01-appb-M000030
 であるため、最終的に式(21)は、
Figure JPOXMLDOC01-appb-M000031
 により計算できる。
　一方、Ｍ＝２のときは
Figure JPOXMLDOC01-appb-M000032
 とすると、
Figure JPOXMLDOC01-appb-M000033
 であるため、最終的に式(21)は、
Figure JPOXMLDOC01-appb-M000034
 により計算できる。
　以上により、(i)ＹＴＹを求め、(ii)ｔｒＹＴＹ，ｄｅｔＹＴＹを求め、(iii)σ１，σ２を求め、(iv)α，βを求め、(v)ｐｒｏｘ｜｜・｜｜＊，ｗ（Ｙ）を以上の式で求めるという手順で加重核型ノルムのＳＶＴを求めることができる。(i), (ii), (v)は積和計算、(iii)は四則計算と根号計算、(iv)は四則計算とクリッピングで算出できる。よって多数の入力行列Ｙ１，…，ＹＬのＳＶＴである
Figure JPOXMLDOC01-appb-M000035
 をＳＩＭＤ演算によって高速、並列に求めることができる。
　図２は、Ｍ≧３の場合にＳＶＴ（Singular Value Thresholding）を算出するアルゴリズムを示す概略図である。このアルゴリズムは疑似的なコードによって記述されている。以下、この図に沿ってアルゴリズムを説明する。
　本アルゴリズムにおいて、入力は、Ｍ×２（Ｍ行２列）の行列Ｙ、およびＳＶＴの重みベクトルｗである。Ｙの１列目、２列目の列ベクトル（Ｍ次元）をそれぞれｙ１，ｙ２と表す。また、ｗ＝［ｗ１，ｗ２］Ｔはｗ２≧ｗ１≧０である。また出力は行列Ｙの重みｗによるＳＶＴである。出力される行列をＺと表す。
　以下では図の左側に付した行番号を参照しながら説明する。第１行において、変数ａ、ｂ、ｃへの代入が行われる。変数ａにはｙ１ Ｔｙ１、変数ｂにはｙ１ Ｔｙ２、変数ｃにはｙ２ Ｔｙ２を代入する。これらａ、ｂ、ｃはいずれもスカラー値である。
　第２行において、変数ｄ、ｅ、ｆへの代入が行われる。変数ｄにはａｃ－ｂ２の値、変数ｅには√ｄの値、変数ｆにはａ＋ｃの値を代入する。変数ｄの値は行列ＹＴＹの行列式であり、変数ｅの値はその平方根、変数ｆの値は行列ＹＴＹのトレースである。
　第３行において、変数ｇ、ｈへの代入が行われる。変数ｇには√（ｆ＋２ｅ）の値、変数ｈには√（ｆ－２ｅ）の値を代入する。変数ｇの値は行列Ｙの特異値の和であり、変数ｈの値は行列Ｙの特異値の差である。
　第４行において、変数σ１，σ２への代入が行われる。変数σ１には（ｇ＋ｈ）／２の値、変数σ２には（ｇ－ｈ）／２の値を代入する。変数σ１の値は行列Ｙの最大特異値であり、変数σ２の値は行列Ｙの最小特異値である。
　第５行において、変数αへの代入が行われる。変数αに、ｇが０でない場合は
Figure JPOXMLDOC01-appb-M000036
 の値を、ｇが０の場合は０を代入する。
　第６行において、変数βへの代入が行われる。変数βに、ｈが０でない場合は
Figure JPOXMLDOC01-appb-M000037
 の値を、ｈが０の場合は０を代入する。
　第７行において、出力である変数Ｚへの代入が行われる。変数Ｚに
Figure JPOXMLDOC01-appb-M000038
 の値を代入する。そして処理を終了してＺを出力する。
　以上、説明したように、ＳＶＤを用いず、少ない計算量でＳＶＴを算出することができる。
　図３は、Ｍ＝２の場合にＳＶＴ（Singular Value Thresholding）を算出するアルゴリズムを示す概略図である。このアルゴリズムは疑似的なコードによって記述されている。以下、この図に沿ってアルゴリズムを説明する。
　本アルゴリズムにおいて、入力は、２×２（２行２列）の行列Ｙ、およびＳＶＴの重みＷである。Ｙの１列目、２列目の列ベクトル（２次元）をそれぞれｙ１，ｙ２と表し、行列Ｙの各要素を、行番号および列番号をこの順で並べたサフィックス（添え字）を用いて、ｙ１，１，ｙ１，２，ｙ２，１，ｙ２，２で表す。また、ｗ＝［ｗ１，ｗ２］Ｔはｗ２≧ｗ１≧０である。また出力は行列Ｙの重みｗによるＳＶＴである。出力される行列をＺと表す。
　以下では図の左側に付した行番号を参照しながら説明する。第１行において、変数ａ、ｃへの代入が行われる。変数ａにはｙ１ Ｔｙ１、変数ｃにはｙ２ Ｔｙ２を代入する。これらａ、ｃはいずれもスカラー値である。
　第２行において、変数ｄ、ｅ、ｆへの代入が行われる。変数ｄにはｙ１，１ｙ２，２－ｙ１，２ｙ２，１の値、変数ｅには｜ｄ｜の値、変数ｆにはａ＋ｃの値を代入する。変数ｄの値は行列Ｙの行列式であり、変数ｅの値はその絶対値で行列ＹＴＹの行列式の平方根、変数ｆの値は行列ＹＴＹのトレースである。
　第３行において、変数ｇ、ｈへの代入が行われる。変数ｇには√（ｆ＋２ｅ）の値、変数ｈには√（ｆ－２ｅ）の値を代入する。変数ｇの値は行列Ｙの特異値の和であり、変数ｈの値は行列Ｙの特異値の差である。
　第４行において、変数σ１，σ２への代入が行われる。変数σ１には（ｇ＋ｈ）／２の値、変数σ２には（ｇ－ｈ）／２の値を代入する。変数σ１の値は行列Ｙの最大特異値であり、変数σ２の値は行列Ｙの最小特異値である。
　第５行において、変数αへの代入が行われる。変数αに、ｇが０でない場合は
Figure JPOXMLDOC01-appb-M000039
 の値を、ｇが０の場合は０を代入する。
　第６行において、変数βへの代入が行われる。変数βに、ｈが０でない場合は
Figure JPOXMLDOC01-appb-M000040
 の値を、ｈが０の場合は０を代入する。
　第７行において、出力である変数Ｚへの代入が行われる。変数Ｚに
Figure JPOXMLDOC01-appb-M000041
 の値を代入する。そして処理を終了してＺを出力する。
　以上、説明したように、ＳＶＤを用いず、少ない計算量でＳＶＴを算出することができる。
　ここで本実施形態のアルゴリズムによるＳＶＴ算出のための計算量について考察する。行列ＹのサイズをＭ×２とする。ＳＶＤを用いるＳＶＴ算出法（従来技術）では、ＳＶＤを求めるために２４Ｍ＋１６０回、閾値処理に２回、行列積を求めるために６Ｍ＋４回の浮動小数点演算が必要である。即ち合計で、３０Ｍ＋１６６回の浮動小数点演算が必要である。
　一方、本実施形態のアルゴリズム１による方法では１２Ｍ＋２６回、アルゴリズム２による方法では３３回の浮動小数点演算が必要である。従って、Ｍ＝２の場合は約８５％、Ｍ＝３の場合は約７６％、Ｍ＝１０の場合は約６９％、それ以上の場合でも６０％以上の浮動小数点演算を削減できることが分かる。
　加えてアルゴリズム１およびアルゴリズム２のアルゴリズムは、ＳＩＭＤなどの並列アーキテクチャにより効率的に並列処理を行うことができる。
［３．行列単純化装置］
　次に、本実施形態の行列単純化装置１について説明する。図４は、本実施形態による行列単純化装置１の概略機能構成を示すブロック図である。図４に示すように、行列単純化装置１は、入力部１１と、ベクトル化部１２と、低ランク近似化部１３と、出力部１４とを含んで構成される。また、各部は、必要に応じてデータを記憶するための記憶部を内部に備える。この記憶部は、半導体メモリーや磁気ハードディスク装置などといった記憶手段を用いて実現される。また、各部の機能は、コンピュータとプログラムとで実現されてもよい。
　行列単純化装置１は、Ｍ行２列（ただし、Ｍ≧２）または２行Ｍ列の行列と２次元ベクトルである重みｗ＝［ｗ１，ｗ２］Ｔを入力され、その行列を低ランク近似化し、低ランク近似化した行列を出力する装置である。
　入力部１１は、Ｍ行２列または２行Ｍ列の行列のデータを外部から取得する。なお、Ｍは、２以上の整数である。この行列の要素は、数値（スカラー値）である。入力部１１は、必要に応じて、入力した行列を転置する。つまり、後段のベクトル化部１２および低ランク近似化部１３がＭ行２列または２行Ｍ列のいずれか一方の形式の行列のみを処理するように構成されているときであって、入力された行列がその形式に合わないとき（つまり、縦と横が逆）に、入力部１１は、入力された行列を転置する。これにより、行列単純化装置１は、Ｍ行２列または２行Ｍ列のいずれの行列をも処理することができるようになる
　以下において、ベクトル化部１２と低ランク近似化部１３とは、Ｍ行２列の行列を処理するものとして説明する。但し、これが２行Ｍ列の行列を処理するものであってもよく、本質的な処理内容は変わらない。
　ベクトル化部１２は、Ｍ行２列の行列（Ｙとする）をベクトル化する。ここでのベクトル化とは、Ｍ行２列の行列Ｙを、２個のＭ次元の列ベクトルｙ１（第１ベクトル）とｙ２（第２ベクトル）とに分割して出力する処理である。つまり、Ｙ＝［ｙ１，ｙ２］である。ベクトル化部１２は、これらのベクトルｙ１，ｙ２を、低ランク近似化部１３に渡す。
　つまり、ベクトル化部１２は、Ｍ行２列の入力行列から各列に対応するＭ次元の第１ベクトルおよび第２ベクトルを抽出する（または２行Ｍ列の入力行列から各行に対応するＭ次元の第１ベクトルおよび第２ベクトルを抽出する）。
　低ランク近似化部１３は、ベクトル化部１２から渡されたベクトルｙ１，ｙ２と重みｗ＝［ｗ１，ｗ２］Ｔに基づき、行列Ｙを低ランク近似化する。言い換えれば、低ランク近似化部１３は、行列Ｙを単純化する。すなわち、低ランク近似化部１３は、Fast Parallel SVT（ＦＰＳＶＴ）で前述した算出方法により、行列Ｙを低ランク近似化する。具体的には、低ランク近似化部１３は、Ｍの値に応じて、前述のアルゴリズム１または２のいずれかを用いて、行列Ｙの低ランク近似化を行う。より具体的には、Ｍ≧３の場合には、低ランク近似化部１３は、アルゴリズム１を用いる。また、Ｍ＝２の場合には、低ランク近似化部１３は、アルゴリズム２を用いる。これにより、低ランク近似化部１３は、低ランク化された行列Ｚを出力する。
　出力部１４は、低ランク近似化部１３によって求められた行列Ｚのデータを外部に出力する。つまり、行列単純化装置１は、入力される行列Ｙを低ランク単純化し、その結果である行列Ｚを出力する。なお、入力部１１が行列の転置を行った場合、出力部１４は、低ランク近似化部１３によって得られた行列Ｚを再び転置してから出力する。これにより、入力される行列のサイズ（行および列の数）と、出力される行列のサイズとを合わせることができる。
　次に、行列単純化装置１が扱う行列およびベクトルのデータの構造について説明する。
　図５Ａ～図５Ｃは、行列単純化装置１が処理する主要なデータの構成を示す概略図である。図５Ａは、入力行列Ｙのデータ構成を示す。ここでは、Ｍ行２列の場合の行列を示している。この図では行番号と列番号を付して示しており、行列Ｙの要素であるｙｉｊ（ｉ＝１，・・・，Ｍ、ｊ＝１，２）が各領域に格納される。なお、２行Ｍ列の行列の場合には、行と列の方向が入れ替わる。図５Ｂは、行列Ｙを基にベクトル化されたベクトルｙ１およびｙ２のデータ構成を示す。ベクトルｙ１およびｙ２は、それぞれ、Ｍ次元の列ベクトルである。行列Ｙ（図５Ａ）の第１列がベクトルｙ１に対応し、第２列がベクトルｙ２に対応する。図５Ｃは、出力行列Ｚのデータ構成を示す。行列Ｚについても、Ｍ行２列の場合の行列を示している。行列Ｚの要素であるzｉｊ（ｉ＝１，・・・，Ｍ、ｊ＝１，２）が各領域に格納される。
　次に、本実施形態の行列単純化装置１の処理の流れについて説明する。図６は、本実施形態による行列単純化装置１の処理の流れを示すフローチャートである。
　まず、入力部１１は、Ｍ行２列（ただし、Ｍ≧２）または２行Ｍ列の行列と２次元ベクトルである重みｗ＝［ｗ１，ｗ２］Ｔを受付ける（ステップＳ１０１）。
　次に、入力部１１は、入力行列が２行Ｍ列の場合（２列でない場合）は、転置してＭ行２列とする（ステップＳ１０２）。
　次に、ベクトル化部１２は、Ｍ行２列の行列Ｙを、２個のＭ次元の列ベクトルｙ１（第１ベクトル）とｙ２（第２ベクトル）とに分割する（ステップＳ１０３）。
　次に、低ランク近似化部１３は、Ｍ≧３か否か判定する（ステップＳ１０４）。
　ステップＳ１０４でＹｅｓ（Ｍ≧３）の場合は、低ランク近似化部１３は、ｙ１ Ｔｙ１、ｙ１ Ｔｙ２、ｙ２ Ｔｙ２を算出して、それぞれ変数ａ、ｂ、ｃへ代入する（ステップＳ１０５）。これらａ、ｂ、ｃはいずれもスカラー値である。
　次に、低ランク近似化部１３は、ａｃ－ｂ２、√ｄ、ａ＋ｃを算出して、それぞれ変数ｄ、ｅ、ｆへ代入する（ステップＳ１０６）。変数ｄの値は行列ＹＴＹの行列式であり、変数ｅの値はその平方根、変数ｆの値は行列ＹＴＹのトレースである。
　次に、低ランク近似化部１３は、√（ｆ＋２ｅ）、√（ｆ－２ｅ）を算出して、それぞれ変数ｇ、ｈへ代入する（ステップＳ１０７）。変数ｇの値は行列Ｙの特異値の和であり、変数ｈの値は行列Ｙの特異値の差である。
　次に、低ランク近似化部１３は、（ｇ＋ｈ）／２、（ｇ－ｈ）／２を算出して、それぞれ変数σ１，σ２へ代入する（ステップＳ１０８）。変数σ１の値は行列Ｙの最大特異値であり、変数σ２の値は行列Ｙの最小特異値である。
　次に、低ランク近似化部１３は、ｇが０でない場合は
Figure JPOXMLDOC01-appb-M000042
 を算出し、変数αへ代入する。一方、低ランク近似化部１３は、ｇが０の場合は変数αへ０を代入する（ステップＳ１０９）。
　次に、低ランク近似化部１３は、ｈが０でない場合は
Figure JPOXMLDOC01-appb-M000043
 を算出し、変数βへ代入する。一方、低ランク近似化部１３は、ｈが０の場合は変数βへ０を代入する（ステップＳ１１０）。
　次に、低ランク近似化部１３は、
Figure JPOXMLDOC01-appb-M000044
 を算出し、変数Ｚへ代入する（ステップＳ１１１）。
　一方、ステップＳ１０４でＮｏ（Ｍ＝２）の場合は、低ランク近似化部１３は、ｙ１ Ｔｙ１、ｙ２ Ｔｙ２を算出して、それぞれ変数ａ、ｃへ代入する（ステップＳ１１２）。これらａ、ｃはいずれもスカラー値である。
　次に、低ランク近似化部１３は、ｙ１，１ｙ２，２－ｙ１，２ｙ２，１、｜ｄ｜、ａ＋ｃを算出して、それぞれ変数ｄ、ｅ、ｆへ代入する（ステップＳ１１３）。変数ｄの値は行列Ｙの行列式であり、変数ｅの値はその絶対値で行列ＹＴＹの行列式の平方根、変数ｆの値は行列ＹＴＹのトレースである。
　次に、低ランク近似化部１３は、√（ｆ＋２ｅ）、√（ｆ－２ｅ）を算出して、それぞれ変数ｇ、ｈへ代入する（ステップＳ１１４）。変数ｇの値は行列Ｙの特異値の和であり、変数ｈの値は行列Ｙの特異値の差である。
　次に、低ランク近似化部１３は、（ｇ＋ｈ）／２、（ｇ－ｈ）／２を算出して、それぞれ変数σ１，σ２へ代入する（ステップＳ１１５）。変数σ１の値は行列Ｙの最大特異値であり、変数σ２の値は行列Ｙの最小特異値である。
　次に、低ランク近似化部１３は、ｇが０でない場合は
Figure JPOXMLDOC01-appb-M000045
 を算出し、変数αへ代入する。一方、低ランク近似化部１３は、ｇが０の場合は変数αへ０を代入する（ステップＳ１１６）。
　次に、低ランク近似化部１３は、ｈが０でない場合は
Figure JPOXMLDOC01-appb-M000046
 を算出し、変数βへ代入する。一方、低ランク近似化部１３は、ｈが０の場合は変数βへ０を代入する（ステップＳ１１７）。
　次に、低ランク近似化部１３は、
Figure JPOXMLDOC01-appb-M000047
 を算出し、変数Ｚへ代入する（ステップＳ１１８）。
　ステップＳ１１１又はステップＳ１１８の後、出力部１４は、変数Ｚ（ステップＳ１０２で入力行列を転置していた場合は、ＺＴ）の値を出力する。
　以上で、図６の行列単純化装置１の処理の流れの説明は終了である。
　上述した各実施形態における行列単純化装置１の機能の全部または一部を、コンピュータで実現するようにしても良い。その場合、この機能を実現するためのプログラムをコンピュータ読み取り可能な記録媒体に記録して、この記録媒体に記録されたプログラムをコンピュータシステムに読み込ませ、実行することによって実現しても良い。なお、ここでいう「コンピュータシステム」とは、ＯＳや周辺機器等のハードウェアを含むものとする。また、「コンピュータ読み取り可能な記録媒体」とは、フレキシブルディスク、光磁気ディスク、ＲＯＭ、ＣＤ－ＲＯＭ、ＤＶＤ－ＲＯＭ、ＵＳＢメモリー等の可搬媒体、コンピュータシステムに内蔵されるハードディスク等の記憶装置のことをいう。さらに「コンピュータ読み取り可能な記録媒体」とは、インターネット等のネットワークや電話回線等の通信回線を介してプログラムを送信する場合の通信線のように、短時間の間、動的にプログラムを保持するもの、その場合のサーバーやクライアントとなるコンピュータシステム内部の揮発性メモリーのように、一定時間プログラムを保持しているものも含んでも良い。また上記プログラムは、前述した機能の一部を実現するためのものであっても良く、さらに前述した機能をコンピュータシステムにすでに記録されているプログラムとの組み合わせで実現できるものであっても良い。
［効果］
　以上、説明した実施形態による行列単純化装置１は、特異値閾値処理（ＳＶＴ）の高速算出アルゴリズムを実行可能である。これにより、多数の小型行列を低ランク正則化する最適化問題の解を高速に求めることが可能となる。
　具体的には、行列単純化装置１は、Ｍ行２列または２行Ｍ列（Ｍ≧２）の行列を対象として、低ランク化処理を、少ない計算量で高速に実行することが可能となる。また、行列単純化装置１は、複数の行列を対象として、低ランク化処理を並列して実行することができる。つまり、行列単純化装置１により、Ｍ行２列または２行Ｍ列の行列に対するＦＰＳＶＴ（Fast Parallel SVT）アルゴリズムが実現される。
　さらに具体的には、核型ノルムを部分空間上のベクトル距離で表現できる発見に基づき、行列単純化装置１により、特異値分解（ＳＶＤ）が不要でデータ並列なＳＶＴ算出法を実現した。また、実際のデータを用いた評価実験の結果、従来手法以上の計算精度を持ちつつ、最大７７倍高速にＳＶＴを算出できることを確認した。
　以上、本発明の実施形態について図面を参照して詳述してきたが、具体的な構成はこの実施形態に限られるものではなく、本発明の要旨を逸脱しない範囲の設計変更等も含まれる。
　本発明は、多数の小型行列が加重核型ノルムで正則化される最適化問題を取り扱う技術に適用できる。またこれにより本発明は、グラフ単純化や偽色除去等をはじめとするデータや物理現象に潜在する低ランク性に基づくデータ解析を行う技術に適用できる。
　１・・・・行列単純化装置
　１１・・・入力部
　１２・・・ベクトル化部
　１３・・・低ランク近似化部
　１４・・・出力部

Claims (8)

  1. 　Ｍ行２列行列Ｙ＝［ｙ１，ｙ２］又は２行Ｍ列行列ＹＴ＝［ｙ１，ｙ２］Ｔを、２つのＭ次元ベクトルｙ１，ｙ２とするベクトル化部と、
    　重みｗ＝［ｗ１，ｗ２］Ｔによる行列Ｙの加重核型ノルムが、像ＩｍＹ上の回転行列をＲとするとき、［ｙ１－Ｒｙ２，ｙ１＋Ｒｙ２］と前記重みｗに基づくグループ尖星準ノルムで表されることに基づいて、前記加重核型ノルムのＳＶＴ（Singular Value Thresholding）を算出することにより、行列Ｙの低ランク近似化を行う低ランク近似化部と、
    　を備える行列単純化装置。
  2. 　前記グループ尖星準ノルムは、
    Figure JPOXMLDOC01-appb-M000001
     とするとき、
    Figure JPOXMLDOC01-appb-M000002
     と表される
    　請求項１に記載の行列単純化装置。
  3. 　前記低ランク近似化部は、前記重みｗによる前記行列Ｙの前記加重核型ノルムのＳＶＴを、前記行列Ｙの最大特異値をσ１、最小特異値をσ２とし、σ１とσ２とｗとに基づく値をα，βとするとき、
    Figure JPOXMLDOC01-appb-M000003
     により算出する
    　請求項１に記載の行列単純化装置。
  4. 　前記低ランク近似化部は、（・）＋をランプ関数とするとき、前記α，βを
    Figure JPOXMLDOC01-appb-M000004
     により算出する
    　請求項３に記載の行列単純化装置。
  5. 　前記低ランク近似化部は、Ｍ≧３のとき、ｄｅｔＹをＹの行列式とするとき、Ｒ［ｙ２，－ｙ１］を
    Figure JPOXMLDOC01-appb-M000005
     により算出する
    　請求項３に記載の行列単純化装置。
  6. 　前記低ランク近似化部は、Ｍ＝２のとき、
    Figure JPOXMLDOC01-appb-M000006
     とし、ｓｇｎ（・）を符号関数、ｄｅｔＹをＹの行列式とするとき、Ｒ［ｙ２，－ｙ１］を
    Figure JPOXMLDOC01-appb-M000007
     により算出する
    　請求項３に記載の行列単純化装置。
  7. 　Ｍ行２列行列Ｙ＝［ｙ１，ｙ２］又は２行Ｍ列行列ＹＴ＝［ｙ１，ｙ２］Ｔを、２つのＭ次元ベクトルｙ１，ｙ２とするベクトル化ステップと、
    　行列Ｙの重みｗ＝［ｗ１，ｗ２］Ｔによる加重核型ノルムが、像ＩｍＹ上の回転行列をＲとするとき、（ｙ１－Ｒｙ２，ｙ１＋Ｒｙ２）とｗに基づくグループ尖星準ノルムで表されることに基づいて、前記加重核型ノルムのＳＶＴ（Singular Value Thresholding）を算出することにより、行列Ｙの低ランク近似化を行う低ランク近似化ステップと、
    　を含む行列単純化方法。
  8. 　コンピュータを、
    　Ｍ行２列行列Ｙ＝［ｙ１，ｙ２］又は２行Ｍ列行列ＹＴ＝［ｙ１，ｙ２］Ｔを、２つのＭ次元ベクトルｙ１，ｙ２とするベクトル化部、
    　行列Ｙの重みｗ＝［ｗ１，ｗ２］Ｔによる加重核型ノルムが、像ＩｍＹ上の回転行列をＲとするとき、（ｙ１－Ｒｙ２，ｙ１＋Ｒｙ２）とｗに基づくグループ尖星準ノルムで表されることに基づいて、前記加重核型ノルムのＳＶＴ（Singular Value Thresholding）を算出することにより、行列Ｙの低ランク近似化を行う低ランク近似化部、
    　として機能させるためのプログラム。
PCT/JP2022/027010 2022-07-07 2022-07-07 行列単純化装置、行列単純化方法、およびプログラム WO2024009469A1 (ja)

Priority Applications (1)

Application Number Priority Date Filing Date Title
PCT/JP2022/027010 WO2024009469A1 (ja) 2022-07-07 2022-07-07 行列単純化装置、行列単純化方法、およびプログラム

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
PCT/JP2022/027010 WO2024009469A1 (ja) 2022-07-07 2022-07-07 行列単純化装置、行列単純化方法、およびプログラム

Publications (1)

Publication Number Publication Date
WO2024009469A1 true WO2024009469A1 (ja) 2024-01-11

Family

ID=89453132

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2022/027010 WO2024009469A1 (ja) 2022-07-07 2022-07-07 行列単純化装置、行列単純化方法、およびプログラム

Country Status (1)

Country Link
WO (1) WO2024009469A1 (ja)

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
ANONYMOUS: "Takayuki Sasaki", 18 April 2023 (2023-04-18), XP093126909, Retrieved from the Internet <URL:https://takasasa.github.io/> *
SASAKI TAKAYUKI, KITAHARA MASAKI, SHIMIZU ATSUSHI: "Mathematical Derivation of Fast Singular Value Thresholding for Low-rank Optimization", FORUM ON INFORMATION TECHNOLOGY 2017 (FIT), PART 1, 1 January 2017 (2017-01-01), pages 5 - 12, XP093126903 *

Similar Documents

Publication Publication Date Title
CN109919183B (zh) 一种基于小样本的图像识别方法、装置、设备及存储介质
CN111295675B (zh) 用于使用内核来处理卷积运算的设备和方法
Nguyen et al. A flexible tensor block coordinate ascent scheme for hypergraph matching
Dhodhi et al. D-ISODATA: A distributed algorithm for unsupervised classification of remotely sensed data on network of workstations
WO2022006919A1 (zh) 基于激活定点拟合的卷积神经网络训练后量化方法及系统
EP3583575B1 (en) Image transformation for machine learning
US20180129930A1 (en) Learning method based on deep learning model having non-consecutive stochastic neuron and knowledge transfer, and system thereof
CN110334761B (zh) 基于正交性约束增量非负矩阵分解的有监督图像识别方法
Sun et al. Multiobjective task scheduling for energy-efficient cloud implementation of hyperspectral image classification
Lee et al. Learning Shared Knowledge for Deep Lifelong Learning using Deconvolutional Networks.
Anderson et al. An effective decoupling method for matrix optimization and its application to the ICA problem
JP6810003B2 (ja) 行列単純化装置、プログラム、および行列単純化方法
Koehl et al. Statistical physics approach to the optimal transport problem
Sukhorukova et al. A generalisation of de la Vallée-Poussin procedure to multivariate approximations
JPWO2016125500A1 (ja) 特徴変換装置、認識装置、特徴変換方法及びコンピュータ読み取り可能記録媒体
WO2024009469A1 (ja) 行列単純化装置、行列単純化方法、およびプログラム
Zerrouk et al. Evolutionary algorithm for optimized CNN architecture search applied to real-time boat detection in aerial images
JP3809062B2 (ja) マルチレベル不完全ブロック分解による前処理を行う処理装置
JP2021005242A (ja) 情報処理装置、情報処理プログラム、及び情報処理方法
Yang et al. Adaptive DNN Surgery for Selfish Inference Acceleration with On-demand Edge Resource
EP3923199A1 (en) Method and system for compressing a neural network
KR20230044830A (ko) 이미지 처리 방법 및 장치
WO2020040007A1 (ja) 学習装置、学習方法及び学習プログラム
Li et al. Hessian Aware Low-Rank Perturbation for Order-Robust Continual Learning
Lin et al. Fast and Sample-Efficient Relevance-Based Multi-Task Representation Learning

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 22950263

Country of ref document: EP

Kind code of ref document: A1