WO2020162189A1 - 非負値行列分解最適化装置、非負値行列分解最適化方法、プログラム - Google Patents

非負値行列分解最適化装置、非負値行列分解最適化方法、プログラム Download PDF

Info

Publication number
WO2020162189A1
WO2020162189A1 PCT/JP2020/002206 JP2020002206W WO2020162189A1 WO 2020162189 A1 WO2020162189 A1 WO 2020162189A1 JP 2020002206 W JP2020002206 W JP 2020002206W WO 2020162189 A1 WO2020162189 A1 WO 2020162189A1
Authority
WO
WIPO (PCT)
Prior art keywords
matrix
optimization
negative
updating
following equation
Prior art date
Legal status (The legal status 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 status listed.)
Ceased
Application number
PCT/JP2020/002206
Other languages
English (en)
French (fr)
Inventor
健太 丹羽
登 原田
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
NTT Inc
Original Assignee
Nippon Telegraph and Telephone Corp
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 Nippon Telegraph and Telephone Corp filed Critical Nippon Telegraph and Telephone Corp
Priority to US17/428,263 priority Critical patent/US20220121728A1/en
Publication of WO2020162189A1 publication Critical patent/WO2020162189A1/ja
Anticipated expiration legal-status Critical
Ceased legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; 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/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; 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
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; 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/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N99/00—Subject matter not provided for in other groups of this subclass

Definitions

  • the present invention relates to a non-negative matrix decomposition technique that represents a non-negative matrix as a product of two non-negative matrices.
  • Non-negative matrix decomposition problem the problem of estimating the non-negative matrix group that potentially constitutes the matrix from the non-negative matrix (hereinafter referred to as the observation matrix) that represents the observed signal or the observed data (non-negative matrix decomposition problem). For example, it is being practically used for a problem of removing noise mixed with a voice spectrogram or an image as an input, and a problem of decomposing and extracting an inherent user characteristic by inputting statistically aggregated data such as a questionnaire. .. In addition, in such application examples, there are increasing situations where reduction of processing time and real-time processing are required.
  • ⁇ represents a superscript, for example, x y ⁇ z represents that y z is a superscript for x.
  • (Underscore) represents a subscript, for example, xy_z represents that y z is a subscript for x.
  • the superscript “ ⁇ ” or “ ⁇ ” such as ⁇ x or ⁇ x for a certain character x should be written directly above “x”, but ⁇ x or ⁇ x Will be described.
  • the non-negative matrix Z is, for example, a matrix representing a noise-removed speech spectrogram or a noise-removed image.
  • the non-negative matrix Z is represented by the product of two non-negative matrices ⁇ A, B ⁇ (where A ⁇ R I ⁇ K and B ⁇ R J ⁇ K ) as in the following equation.
  • K is a sufficiently small number, and K is preferably a number smaller than I,J.
  • T represents transposition.
  • the matrix ⁇ A, B ⁇ may also be called a factorized element matrix.
  • the cost function J(A, B):R I ⁇ K ⁇ R J ⁇ K ⁇ R ⁇ ⁇ - ⁇ , + ⁇ is the normalization term H A ,H for the loss term G and the matrix ⁇ A, B ⁇ . Composed of the sum of B.
  • the loss term G is to measure the errors of the two matrices (in the formula (1-2), the matrix Y and the matrix AB T).
  • the normalization terms H A and H B guarantee the non-negativeness of the matrix ⁇ A, B ⁇ .
  • the problem to be solved is to estimate the matrix ⁇ A, B ⁇ that minimizes the cost function J, which can be formulated by the following equation.
  • Expression (1-3) represents an operation of estimating a nonnegative matrix ⁇ A, B ⁇ that is the lower bound (infimum) of the cost function J.
  • loss term G any closed true convex function (hereinafter simply referred to as a convex function) can be used.
  • the loss term G is expressed by Bregman divergence (which is a convex function in a narrow sense), it is as follows.
  • the loss term G ⁇ ( ⁇ ) can be used as the formula (1-6 ) ⁇ divergence.
  • ⁇ divergence is one of the loss terms often used in nonnegative matrix decomposition.
  • the normalization term it is necessary to guarantee that at least the matrix ⁇ A, B ⁇ treated as a variable is nonnegative.
  • the normalization term can be defined using a non-negative barrier function that outputs 0 when all elements of the matrix are 0 or more and outputs + ⁇ in all other cases.
  • the matrix ⁇ A, B ⁇ has sparseness as the probabilistic property of the matrix ⁇ A, B ⁇ , it can be used as the normalization term of the following equation using the L1 norm.
  • ⁇ A and ⁇ B are the (positive) constants given in advance.
  • Non-Patent Document 1 an optimization algorithm based on the multiplicative update rule (Non-Patent Document 1) and an optimization algorithm based on the Alternating direction method (of Non-Patent Document 2) will be described. ..
  • the multiplicative updating method can be expressed as a kind of the forward gradient method, and is configured by subtracting an amount obtained by multiplying the gradient of the cost function J by the step size ⁇ from the current element.
  • Equations (1-10) and (1-11) are as follows.
  • [z] + max ⁇ z, 0 ⁇ , and ⁇ (>0) is a constant.
  • the observation matrix Y is input, and the matrix ⁇ A, B ⁇ is updated a predetermined number of times by the equations (1-14) and (1-15) to obtain the matrix ⁇ A, B ⁇ .
  • the matrix E representing noise can also be obtained by using the equation (1-1).
  • the non-negative matrix decomposition problem is formulated as an optimization problem of the following formula.
  • Equation (1-18) is the lower bound (infimum) of the Lagrangian function L, and after estimating the matrix Z, A, B, the upper bound (supremum) of the Lagrangian function L for the estimated matrix Z, A, B Represents the operation of estimating the dual variable U.
  • the ADMM-based optimization algorithm for solving the problem of equation (1-18) is as shown in Figs.
  • the optimization algorithm of FIG. 1 and the optimization algorithm of FIG. 2 are equivalent.
  • Variables ⁇ X, ⁇ V ⁇ appearing in the optimization algorithms of FIGS. 1 and 2 are variables obtained by nonlinearly transforming the dual variable U. By using this variable ⁇ X, ⁇ V ⁇ , the formula in the algorithm is simplified.
  • ⁇ (>0) is a constant representing the step size.
  • the optimization algorithms of FIGS. 1 and 2 each take the observation matrix Y as an input and update the matrix ⁇ A, B ⁇ a predetermined number of times to obtain the matrix ⁇ A, B ⁇ .
  • the optimization algorithm based on the multiplicative updating method has the merit that the updating rule is very simple and easy to implement. Further, there is an advantage that the convergence speed is very fast in the initial stage of updating and the calculation amount per step is low. However, it has been reported that the error convergence rate deteriorates as the number of updates increases and approaches a fixed point (that is, a solution that minimizes the cost function) (see Non-Patent Document 2). That is, the optimization algorithm based on the multiplicative updating method has a problem that it is difficult to approach a solution with high estimation accuracy and the convergence is slow.
  • Non-Patent Document 2 it has been experimentally shown that the optimization algorithm based on ADMM tends not to deteriorate the error convergence rate even if the number of updates increases (see Non-Patent Document 2).
  • the step size ( ⁇ in the algorithm of FIGS. 1 and 2) required for updating must be determined in advance, and that the effect cannot be obtained unless the step size is set appropriately.
  • the ADMM-based optimization algorithm has a problem in that convergence may not be stable because it may be difficult to approach a solution with high estimation accuracy such as when the step size is not set appropriately.
  • an object of the present invention is to provide a non-negative matrix decomposing technique that converges rapidly and stably.
  • Y is a non-negative observation matrix representing an observation signal or observation data
  • a non-negative matrix factorization optimizing device that includes an optimizing unit that optimizes a non-negative matrix ⁇ A, B ⁇ that is a factorization
  • a non-negative matrix decomposition that represents a non-negative matrix as a product of two non-negative matrices can be obtained quickly and stably.
  • FIG. 6 is a flowchart showing the operation of the non-negative matrix decomposition optimization apparatus 100.
  • 3 is a block diagram showing a configuration of an optimization unit 110.
  • FIG. 6 is a flowchart showing the operation of the optimization unit 110.
  • ⁇ (Caret) represents a superscript.
  • x y ⁇ z means that y z is a subscript to x
  • x y ⁇ z means that y z is a subscript to x
  • _ (underscore) represents a subscript.
  • x y_z means that y z is a subscript to x
  • x y_z means that y z is a subscript to x.
  • ADMM uses Douglas-Rachford type monotonic operator decomposition. Further, as described above, the ADMM-based optimization algorithm needs to appropriately set predetermined parameters such as a step size in advance.
  • the optimization algorithm in the embodiment of the present invention is configured as an optimization algorithm that converges fast and stably based on the optimization algorithm based on ADMM.
  • the optimization algorithm in the embodiment of the present invention is configured as an optimization algorithm that (1) does not require parameter adjustment such as step size and (2) can also use Peaceman-Rachford type monotonic operator decomposition. ..
  • G * and H * are convex conjugate functions (closed true convex functions) represented by the following equation.
  • T i is a monotonic operator that inputs the I ⁇ J matrix and outputs the I ⁇ J matrix.
  • ⁇ represents composition of operators. Therefore, in equation (2-6), ⁇ D -1 and T i are combined.
  • X ⁇ R I ⁇ J is an auxiliary variable of the dual variable U
  • ⁇ is an averaging coefficient of the averaging operator.
  • ⁇ can be any number that satisfies ⁇ ⁇ (0, 1), but often uses 1/2.
  • the BDR-type monotonic operator decomposition is an application of the averaging operator to the BPR-type monotonic operator decomposition. Is expected to converge rapidly.
  • Formula (2-9) to Formula (2-12) is derived from Formula (2-7), BPR type variable update rule, Formula (2-9) to Formula (2-11) and Formula (2-13) Is the BDR type variable update rule derived from Eq. (2-8).
  • equations (2-9) to (2-13) include the sub-differential operator T i of the convex conjugate function, the equations are further changed.
  • the equation (2-9) is transformed so that the main variable and the dual variable are alternately updated.
  • the problem related to the convex conjugate function G * (see Eq. (2-2)) included in Eq. (2-9) is defined below.
  • variable ⁇ Z, U ⁇ is updated so that the subordinate derivative of the Lagrangian function Q contains the zero vector.
  • variable U is also updated using the updated variable Z.
  • equations (2-9) and (2-10) are expressed as follows.
  • Bregman divergence B D ⁇ + is a strictly convex function that modifies the metric of the variable space, and is given by
  • equations (2-11) to (2-13) are expressed as follows.
  • Expressions (2-14) to (2-19) are expressions that form the variable update rule of the optimization algorithm of the embodiment of the present invention.
  • 3 and 4 show the B-P-R type optimization algorithm and the B-D-R type optimization algorithm.
  • z i,j represents the (i,j) element of the matrix Z.
  • ⁇ (>0) is a sufficiently small number and is used to guarantee that the function D is a strictly convex function.
  • z ij old represents the (i,j) element of the current matrix Z.
  • h ij be the (i,j) element of the matrix H of the following equation.
  • Equation (2-23) shows that (I ⁇ J) matrix elements are represented by (I+J) vector elements, and in addition to speeding up variable updating, it also reduces memory capacity. Contribute.
  • the vector ⁇ h A , h B ⁇ can be determined, for example, by using the average value of the quadratic gradient of the predetermined number of functions G.
  • Y be a non-negative observation matrix representing observation signals or observation data
  • Z be a non-negative matrix forming the observation matrix Y that satisfies equation (1-1).
  • FIG. 5 is a block diagram showing the configuration of the non-negative matrix factorization optimization device 100.
  • FIG. 6 is a flowchart showing the operation of the non-negative matrix factorization optimization device 100.
  • the non-negative matrix decomposition optimization device 100 includes an optimization unit 110 and a recording unit 190.
  • the recording unit 190 is a component that appropriately records information necessary for the processing of the non-negative matrix factorization optimization device 100.
  • the recording unit 190 records, for example, an observation matrix Y that is an input and a matrix ⁇ A, B ⁇ that is an optimization target.
  • the cost function J(A, B): R I ⁇ K ⁇ R J ⁇ K ⁇ R ⁇ - ⁇ , + ⁇ is defined by the following equation.
  • H A is a normalization term for matrix A
  • H B is a normalization term for matrix B.
  • FIG. 7 is a block diagram showing the configuration of the optimization unit 110.
  • FIG. 8 is a flowchart showing the operation of the optimizing unit 110.
  • the optimization unit 110 includes an initialization unit 111, a matrix Z updating unit 112, a first dual variable updating unit 113, a matrix A updating unit 114, a matrix B updating unit 115, and a second dual pair. It includes a variable updating unit 116, a counter updating unit 117, and a termination condition determining unit 118.
  • the matrix Z updating unit 112 updates the matrix Z by the following formula.
  • H Z is the normalization term for the matrix Z
  • B D ⁇ + is the Bregman divergence defined using the function D + .
  • C is a matrix with c ij as the (i, j) element
  • h ij ⁇ ij 2 G(c ij old )( c ij old is the (i, j) element of the current matrix C
  • ⁇ (> (0) is given by a predetermined constant.
  • Z or ⁇ X is used when calculating the function D + .
  • This function D + is also used in the matrix A updating unit 114 and the matrix B updating unit 115. In this case, when calculating the function D +, AB T and - ⁇ V is used.
  • the first dual variable updating unit 113 updates the dual variable ⁇ V by the following equation.
  • the matrix A updating unit 114 updates the matrix A by the following formula.
  • the matrix B updating unit 115 updates the matrix B by the following formula.
  • the second dual variable updating unit 116 updates the dual variable ⁇ X by the following equation.
  • the counter updating unit 117 increments the counter t by 1. Specifically, t ⁇ t+1.
  • the termination condition determination unit 118 determines that the counter t has reached the predetermined number of updates T (T is an integer of 1 or more, for example, 100,000) (that is, t>T, and the termination condition is satisfied). Case) outputs the matrix ⁇ A, B ⁇ at that time and ends the process. Otherwise, the process returns to S112. That is, the optimization unit 110 repeats the processing of S112 to S118.
  • S114 and S115 (that is, the updating process of the matrix A and the updating process of the matrix B) may be performed in the order of S114 and S115, or may be performed in the order of S115 and S114. Further, S114 and S115 may be executed in parallel.
  • the second dual variable updating unit 116 updates the dual variable ⁇ X according to the BPR type updating rule, it may be updated according to the BDR type updating rule. In this case, in S116, the second dual variable updating unit 116 updates the dual variable ⁇ X by the following equation.
  • the optimizing unit 110 optimizes the matrix ⁇ A, B ⁇ based on the Bregman monotonic operator decomposition.
  • a non-negative matrix decomposition that represents a non-negative matrix as a product of two non-negative matrices can be obtained quickly and stably.
  • the function D + is adaptively modified according to the convex function G
  • the metric of the variable space can be adaptively modified, which is necessary in the conventional ADMM-based optimization algorithm (for example, step No need to adjust parameters such as size.
  • the device of the present invention is, for example, as a single hardware entity, an input unit to which a keyboard or the like can be connected, an output unit to which a liquid crystal display or the like can be connected, and a communication device (for example, a communication cable) capable of communicating outside the hardware entity.
  • a communication device for example, a communication cable
  • Connectable communication unit CPU (Central Processing Unit, cache memory and registers may be provided), memory RAM and ROM, hard disk external storage device and their input unit, output unit, communication unit , A CPU, a RAM, a ROM, and a bus connected so that data can be exchanged among external storage devices.
  • the hardware entity may be provided with a device (drive) capable of reading and writing a recording medium such as a CD-ROM.
  • a device having such hardware resources there is a general-purpose computer or the like.
  • the external storage device of the hardware entity stores a program necessary to realize the above-described functions and data necessary for the processing of this program (not limited to the external storage device, for example, the program is read). It may be stored in a ROM that is a dedicated storage device). In addition, data and the like obtained by the processing of these programs are appropriately stored in the RAM, the external storage device, or the like.
  • each program stored in an external storage device or ROM, etc.
  • data necessary for the processing of each program are read into the memory as necessary, and interpreted and executed/processed by the CPU as appropriate. ..
  • the CPU realizes a predetermined function (each constituent element represented by the above,... Unit,... Means, etc.).
  • the present invention is not limited to the above-mentioned embodiments, and can be modified as appropriate without departing from the spirit of the present invention. Further, the processes described in the above-described embodiments are not only executed in time series in the order described, but may be executed in parallel or individually according to the processing capability of the device that executes the processes or as necessary. ..
  • the processing function of the hardware entity (the device of the present invention) described in the above embodiment is realized by a computer
  • the processing content of the function that the hardware entity should have is described by the program.
  • the processing functions of the hardware entity are realized on the computer.
  • the program describing this processing content can be recorded in a computer-readable recording medium.
  • the computer-readable recording medium may be, for example, a magnetic recording device, an optical disc, a magneto-optical recording medium, a semiconductor memory, or the like.
  • a magnetic recording device a hard disk device, a flexible disk, a magnetic tape, etc.
  • an optical disc a DVD (Digital Versatile Disc), a DVD-RAM (Random Access Memory), a CD-ROM (Compact Disc Read Only). Memory), CD-R (Recordable)/RW (ReWritable), etc. as magneto-optical recording medium, MO (Magneto-Optical disc) etc. as semiconductor memory, EEP-ROM (Electronically Erasable and Programmable-Read Only Memory) etc. Can be used.
  • Distribution of this program is carried out by selling, transferring, or lending a portable recording medium such as a DVD or a CD-ROM in which the program is recorded. Further, the program may be stored in a storage device of a server computer and transferred from the server computer to another computer via a network to distribute the program.
  • a computer that executes such a program first stores, for example, the program recorded in a portable recording medium or the program transferred from the server computer in its own storage device. Then, when executing the process, this computer reads the program stored in its own storage device and executes the process according to the read program.
  • a computer may directly read the program from a portable recording medium and execute processing according to the program, and the program is transferred from the server computer to this computer. Each time, the processing according to the received program may be sequentially executed.
  • ASP Application Service Provider
  • the program in this embodiment includes information that is used for processing by an electronic computer and that is equivalent to the program (data that is not a direct command to a computer, but has the property of defining computer processing).
  • the hardware entity is configured by executing a predetermined program on the computer, but at least a part of the processing contents may be realized by hardware.

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Computing Systems (AREA)
  • Operations Research (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Evolutionary Biology (AREA)
  • Probability & Statistics with Applications (AREA)
  • Complex Calculations (AREA)

Abstract

高速で安定的に収束する非負値行列分解技術を提供する。Yを観測信号または観測データを表す非負の観測行列、Zを観測行列Yを構成する非負の行列とし、観測行列Yを入力として、Z=ABTを満たす行列Zの分解である非負の行列{A, B}を最適化する最適化部を含む非負値行列分解最適化装置であって、最適化部が最適化に用いるコスト関数Jは、式J(A, B)=G(Y||ABT)+HA(A)+HB(B)（ただし、Gは損失項、HAは行列Aに対する正規化項、HBは行列Bに対する正規化項）で定義されるものであり、最適化部は、ブレグマン単調作用素分解に基づいて行列{A, B}を最適化する。

Description

非負値行列分解最適化装置、非負値行列分解最適化方法、プログラム
　本発明は、非負値行列を２つの非負値行列の積として表す非負値行列分解技術に関する。
　観測信号や観測データを表す非負値行列（以下、観測行列という）から当該行列を潜在的に構成する非負値行列群を推定する問題（非負値行列分解問題）が汎用的に利用されている。例えば、音声スペクトログラムまたは画像を入力として混在するノイズを除去する問題や、アンケート等の統計集約データを入力として内在するユーザー特性を分解して抽出する問題（トピックモデル）で実用的に利用されつつある。そして、こうした応用事例において処理時間の軽減やリアルタイム性が求められる場面も増えてきている。
　ここでは、まず、上記非負値行列分解問題の定式化を行う。なお、明細書の記載表記の制約上、次のような表記に関する約束を導入する。^（キャレット）は上付き添字を表すこととし、例えば、xy^zはyzがxに対する上付き添字であることを表す。（アンダースコア）は下付き添字を表すこととし、例えば、xy_zはyzがxに対する上付き添字であることを表す。また、_ある文字xに対する^xや~xのような上付き添え字の”^”や”~”は、本来”x”の真上に記載されるべきであるが、^xや~xと記載することとする。
　（非負の）観測行列Y∈RI×Jは、非負値行列分解可能な行列Z∈RI×Jとノイズを表す行列E∈RI×Jの和で構成されるとする。ここで、非負値行列Zは、例えば、ノイズを除去した音声スペクトログラムまたはノイズを除去した画像を表す行列となる。また、非負値行列Zは次式のように2つの非負値行列{A, B}（ただし、A∈RI×K, B∈RJ×K）の積で表されることとする。
Figure JPOXMLDOC01-appb-M000026
ここで、Kは十分小さな数であり、KはI, Jより小さい数であるのが好ましい。また、Tは転置を表す。なお、行列{A, B}を分解要素行列と呼ぶこともある。
　また、コスト関数J(A, B):RI×K×RJ×K→R∪{-∞, +∞}は、損失項Gと行列{A, B}に対する正規化項HA, HBの和で構成される。
Figure JPOXMLDOC01-appb-M000027
ここで、損失項Gは2つ行列（式(1-2)では、行列Yと行列ABT）の誤差を測定するものである。また、正規化項HA, HBは、行列{A, B}の非負性を保証するものである。
　以上より、解くべき問題は、コスト関数Jを最小化するような行列{A, B}を推定することであり、次式により定式化できる。
Figure JPOXMLDOC01-appb-M000028
ここで、式(1-3)はコスト関数Jの下界(infimum)となる非負値行列{A, B}を推定する操作を表す。
　以下、損失項と正規化項について詳細に説明する。まず、損失項について説明する。損失項Gには、任意の閉真凸関数（以後、単に凸関数と表す）を利用できる。例えば、（狭義の凸関数である）ブレグマンダイバージェンスによって損失項Gを表す場合、以下のようになる。
Figure JPOXMLDOC01-appb-M000029
　例えば、式(1-5)で表される（連続で微分可能かつ強凸である）狭義凸関数Φ(β)を用いることにすると、損失項GΦ^(β)として式(1-6)のβダイバージェンスを得る。βダイバージェンスは非負行列分解でよく利用される損失項の一つである。
Figure JPOXMLDOC01-appb-M000030
なお、このβダイバージェンスは、β=1のときフロベニウスノルム、β=0のとき一般化KLダイバージェンス、β=-1のとき板倉－斎藤距離と一致する。
　次に、正規化項について説明する。正規化項に対しては、少なくとも変数として扱う行列{A, B}が非負であることを保証する必要がある。このために、例えば、行列の全要素が0以上である場合に0を出力、その他の場合には+∞を出力する非負バリア関数を用いて正規化項を定義することができる。
Figure JPOXMLDOC01-appb-M000031
ただし、
Figure JPOXMLDOC01-appb-M000032
(0≦i≦I, 0≦j≦J, 0≦k≦K)である。
　さらに、行列{A, B}の確率的な性質として、行列{A, B}がスパース性を持つことを仮定とすると、L1ノルムを用いた次式の正規化項とすることもできる。
Figure JPOXMLDOC01-appb-M000033
ただし、μA, μBは事前に与えられた（正の）定数である。
　続いて、上記非負値行列分解問題を解くための、非負値行列分解の最適化アルゴリズムについて説明する。ここでは、乗法更新法(multiplicative update rule)（非特許文献１）に基づく最適化アルゴリズムと交互方向乗数法(Alternating direction method of multipliers, ADMM)（非特許文献２）に基づく最適化アルゴリズムについて説明する。
　まず、乗法更新法に基づく最適化アルゴリズムについて説明する。乗法更新法は、前向き勾配法の一種として表現できるものであり、現在の要素からコスト関数Jの勾配にステップサイズηを掛け合わせた量を減算するという処理で構成される。
Figure JPOXMLDOC01-appb-M000034
　例えば、式(1-6)のβダイバージェンスを損失項として用いる場合、コスト関数Jの勾配は次式で与えられる。
Figure JPOXMLDOC01-appb-M000035
　ここで、ステップサイズηik, ηjk,を次式で定義すると、
Figure JPOXMLDOC01-appb-M000036
式(1-10)及び式(1-11)は以下のようになる。
Figure JPOXMLDOC01-appb-M000037
ここで、[z]+=max{z, 0}であり、ε(>0)は定数である。
　したがって、乗法更新法に基づく最適化アルゴリズムでは、観測行列Yを入力とし、式(1-14)及び式(1-15)により行列{A, B}を所定の回数更新することにより、行列{A, B}を求める。なお、式(1-1)を用いることにより、ノイズを表す行列Eも求めることができる。
　次に、ADMMに基づく最適化アルゴリズムについて説明する。なお、以下説明する方法ではコスト関数の形式が非特許文献２のそれとは異なるため、アルゴリズムの形式に差異があるものとなるが、本質的には等価なものである。
　このアルゴリズムにおいて重要なアイデアは、式(1-2)の損失項に含まれる行列積ABTを補助変数Zとして置き換え、制約付凸最小化問題として扱う点にある。
Figure JPOXMLDOC01-appb-M000038
　式(1-16)に関連するラグランジュ関数Lを次式で定義する。
Figure JPOXMLDOC01-appb-M000039
ただし、U∈RI×Jは双対変数である。また、<U, Z>=tr(UTZ)である。
　ここで、（主問題を解く代わりに）双対問題を解くことにすると、非負値行列分解問題は次式の最適化問題として定式化される。
Figure JPOXMLDOC01-appb-M000040
ここで、式(1-18)はラグランジュ関数Lの下界(infimum)となる行列Z, A, Bを推定した後、推定した行列Z, A, Bに対するラグランジュ関数Lの上界(supremum)となる双対変数Uを推定する操作を表す。
　式(1-18)の問題を解く、ADMMに基づく最適化アルゴリズムは図１、図２のようになる。図１の最適化アルゴリズムと図２の最適化アルゴリズムは等価なものである。また、図１・図２の最適化アルゴリズムに現れる変数{~X, ~V}は、双対変数Uを非線形変換して得られる変数である。この変数{~X, ~V}を用いることにより、アルゴリズム中の数式が単純化される。κ(>0)はステップサイズを表す定数である。図１・図２の最適化アルゴリズムは、いずれも、観測行列Yを入力とし、行列{A, B}を所定の回数更新することにより、行列{A, B}を求める。
D. D. Lee and H. S. Seung, "Learning the parts of objects by non-negative matrix factorization", Nature, Vol.401, No.6755, pp.788-791, 1999. D. L. Sun and C. Fevotte, "Alternating direction method of multipliers for non-negative matrix factorization with the beta-divergence", in IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP) 2014. IEEE, pp.6201-6205, 2014.
　乗法更新法に基づく最適化アルゴリズムは、更新則が非常に単純で、実装が容易であるというメリットがある。また、収束速度も更新の初期段階ではとても速くかつ１ステップあたりの計算量も低いというメリットもある。しかし、更新回数が増えて不動点（すなわち、コスト関数を最小化する解）に近づくに従い、誤差の収束率が悪化することが報告されている（非特許文献２参照）。つまり、乗法更新法に基づく最適化アルゴリズムには、推定精度の高い解に近づきにくく、収束が遅いという問題がある。
　一方、ADMMに基づく最適化アルゴリズムは、更新回数が増加しても、誤差の収束率が悪化しにくい傾向が実験的に示されている（非特許文献２参照）。しかし、更新に必要となるステップサイズ（図１・図２のアルゴリズムにおけるκ）を事前に決める必要があり、適切にステップサイズを設定しないと、その効果が得られないことが報告されている。つまり、ADMMに基づく最適化アルゴリズムには、ステップサイズが適切に設定されない場合など推定精度の高い解に近づきにくい場合があり、収束が安定しないという問題がある。
　そこで本発明では、高速で安定的に収束する非負値行列分解技術を提供することを目的とする。
　本発明の一態様は、Yを観測信号または観測データを表す非負の観測行列、Zを観測行列Yを構成する非負の行列とし、観測行列Yを入力として、Z=ABTを満たす行列Zの分解である非負の行列{A, B}を最適化する最適化部を含む非負値行列分解最適化装置であって、前記最適化部が最適化に用いるコスト関数Jは、式J(A, B)=G(Y||ABT)+HA(A)+HB(B)（ただし、Gは損失項、HAは行列Aに対する正規化項、HBは行列Bに対する正規化項）で定義されるものであり、前記最適化部は、ブレグマン単調作用素分解に基づいて行列{A, B}を最適化する。
　本発明によれば、非負値行列を２つの非負値行列の積として表す非負値行列分解を高速かつ安定的に得ることができる。
ADMMを用いた最適化アルゴリズムの一例を示す図である。 ADMMを用いた最適化アルゴリズムの一例を示す図である。 本発明の実施形態における最適化アルゴリズム（B-P-R型）の一例を示す図である。 本発明の実施形態における最適化アルゴリズム（B-D-R型）の一例を示す図である。 非負値行列分解最適化装置１００の構成を示すブロック図である。 非負値行列分解最適化装置１００の動作を示すフローチャートである。 最適化部１１０の構成を示すブロック図である。 最適化部１１０の動作を示すフローチャートである。
　以下、本発明の実施の形態について、詳細に説明する。なお、同じ機能を有する構成部には同じ番号を付し、重複説明を省略する。
　各実施形態の説明に先立って、この明細書における表記方法について説明する。
　^（キャレット）は上付き添字を表す。例えば、xy^zはyzがxに対する上付き添字であり、xy^zはyzがxに対する下付き添字であることを表す。また、_（アンダースコア）は下付き添字を表す。例えば、xy_zはyzがxに対する上付き添字であり、xy_zはyzがxに対する下付き添字であることを表す。
　ある文字xに対する^xや~xのような上付き添え字の”^”や”~”は、本来”x”の真上に記載されるべきであるが、明細書の記載表記の制約上、^xや~xと記載しているものである。
＜技術的背景＞
　ADMMではDouglas-Rachford型単調作用素分解を用いている。また、先述したように、ADMMに基づく最適化アルゴリズムはステップサイズなど所定のパラメータを事前に適切に設定する必要がある。
　上記２点に着目し、本発明の実施形態における最適化アルゴリズムを、ADMMに基づく最適化アルゴリズムをベースとした高速かつ安定的に収束する最適化アルゴリズムとして構成する。具体的には、本発明の実施形態における最適化アルゴリズムを、（１）ステップサイズなどのパラメータ調整を不要とし、（２）Peaceman-Rachford型単調作用素分解も利用可能とした最適化アルゴリズムとして構成する。
　まず、本発明の実施形態における最適化アルゴリズムの背景となる理論について説明する。式(1-18)で表される双対問題は、以下のように式変形することにより、凸共役関数（閉真凸関数）の和の最小化問題として扱える。
Figure JPOXMLDOC01-appb-M000041
　ここで、G*, H*は次式で表される凸共役関数（閉真凸関数）である。
Figure JPOXMLDOC01-appb-M000042
　式(2-1)の不動点は、次式のように凸共役関数{G*, H*}の劣微分の和が零行列を含むときに得られる。
Figure JPOXMLDOC01-appb-M000043
ここで、Ti:RI×J→RI×J(i=1, 2)は凸共役関数{G*, H*}の劣微分作用素を表す。すなわち、
Figure JPOXMLDOC01-appb-M000044
となり、Tiは、I×J行列を入力とし、I×J行列を出力する単調作用素となる。
　パラメータ調整を不要とするためのメインとなるアイデアは、凸関数に適応して変数空間の計量を自動的に修正することにある。ここでは、そのための作用素∇D-1を導入する。ただし、関数D:RI×J→RI×Jは（連続で微分可能かつ強凸である）狭義凸関数であり、その微分作用素∇D:RI×J→RI×Jとその逆作用素∇D-1:RI×J→RI×Jは次式を満たすものとする。
Figure JPOXMLDOC01-appb-M000045
　凸関数に適応して変数空間の計量を自動的に修正する方法については後ほど詳述する。
　関数Dが式(2-5)を満たすとき、不動点に関する条件である式(2-4)の両辺に∇D-1を適用することで以下を得る。
Figure JPOXMLDOC01-appb-M000046
ここで、○は作用素の合成を表す。したがって、式(2-6)では、∇D-1とTiを合成している。
　参考非特許文献１に従って式(2-6)を変形すると、式(2-7)のBregman Peaceman-Rachford (B-P-R)型単調作用素分解と式(2-8)のBregman Douglas-Rachford (B-D-R)型単調作用素分解の2種類の単調作用素分解が得られる。
Figure JPOXMLDOC01-appb-M000047
ここで、X∈RI×Jは双対変数Uの補助変数、ξは平均化作用素の平均化係数である。なお、ξにはξ∈(0, 1)を満たす任意の数を用いることができるが、1/2を用いることが多い。
（参考非特許文献１：K. Niwa and W. B. Kleijn, “Bregman monotone operator splitting”, arXiv preprint arXiv:1807.04871, 2018.）
　式(2-7)、式(2-8)をみるとわかるように、B-D-R型単調作用素分解はB-P-R型単調作用素分解に平均化作用素を適用したものであるため、B-P-R型単調作用素分解の方が高速に収束することが予測される。
　式(2-7)及び式(2-8)から再帰的な変数更新則を導出すると、式(2-9)～式(2-13)のような更新処理が導出される。
Figure JPOXMLDOC01-appb-M000048
　式(2-9)～式(2-12)が式(2-7)から導出されるB-P-R型変数更新則、式(2-9)～式(2-11)と式(2-13)が式(2-8)から導出されるB-D-R型変数更新則である。
　式(2-9)～式(2-13)には凸共役関数の劣微分作用素Tiなどが含まれているため、更なる式変更を行う。まず、式(2-9)を主変数と双対変数を交互に更新するように式変形を行う。式(2-9)に含まれる凸共役関数G*（式(2-2)参照）に関連する問題を以下で定義する。
Figure JPOXMLDOC01-appb-M000049
　上記問題のラグランジュ関数Qを次式で定義する。
Figure JPOXMLDOC01-appb-M000050
　ラグランジュ関数Qの劣微分がゼロベクトルを含むように変数{Z, U}は更新される。
Figure JPOXMLDOC01-appb-M000051
であることから、以下を得る。
Figure JPOXMLDOC01-appb-M000052
　この変数Zと変数Uの関係を利用すると、式(2-9)のレゾルヴェント作用素R1の入出力変数{U, X}に対して、以下が成り立つ。
Figure JPOXMLDOC01-appb-M000053
　この式を変形することで、以下を得る。
Figure JPOXMLDOC01-appb-M000054
　これは以下の式を解くことと等価である。
Figure JPOXMLDOC01-appb-M000055
　更新された変数Zを用いて、変数Uも更新される。
Figure JPOXMLDOC01-appb-M000056
　式(2-10)と上記式を組み合わせ、変数Uを除去することで、双対補助変数Vの更新則が得られる。
Figure JPOXMLDOC01-appb-M000057
　式を簡略化するために、以下の双対補助変数{~V, ~X}を導入する。
Figure JPOXMLDOC01-appb-M000058
　双対補助変数{~V, ~X}を用いると、式(2-9)及び式(2-10)は以下のように表される。
Figure JPOXMLDOC01-appb-M000059
ここで、ブレグマンダイバージェンスBD^+は変数空間の計量を修正する狭義凸関数であり、以下で与えられる。
Figure JPOXMLDOC01-appb-M000060
　同様に、式(2-11)～式(2-13)は以下のように表される。
Figure JPOXMLDOC01-appb-M000061
　式(2-14)～式(2-19)が本発明の実施形態の最適化アルゴリズムの変数更新則を構成する式である。図３、図４に、B-P-R型の最適化アルゴリズム、B-D-R型の最適化アルゴリズムを示す。
　最後に、パラメータ調整をしなくても高速に収束することを実現するための狭義凸関数D+を設計する方法について説明する。具体的には、凸関数Gに合わせて関数D+を適応的に修正する方法を説明する。ここでは、一つの実装方法として、次式の一般化二次形式を用いる。なお、計算量低減のために、式(1-4)のように行列要素間の相関を無視したコストを想定している。
Figure JPOXMLDOC01-appb-M000062
ここで、zi,jは行列Zの(i,j)要素を表す。ε(>0)は十分に小さい数であり、関数Dが狭義凸関数であることを保証するために用いる。
　また、∇D+(0)=0を満たすことから、hijは、次式のように、関数Gの二次勾配に従うように設計するのが効率的である。
Figure JPOXMLDOC01-appb-M000063
ただし、zij oldは現在の行列Zの(i,j)要素を表す。
　このとき、関数D+と微分作用素が逆となる関数Dは、次式により与えられる。
Figure JPOXMLDOC01-appb-M000064
　さらに計算量を削減し、高速化するためには、hijが次式の行列Hの(i,j)要素となるようにするとよい。
Figure JPOXMLDOC01-appb-M000065
ここで、hA∈RI, hB∈RJである。
　式(2-23)は（I×J）個の行列要素が（I+J）個のベクトル要素で表されることを示しており、変数更新の高速化に加えて、メモリ量抑制にも寄与する。
　なお、ベクトル{hA, hB}は、例えば、所定の数の関数Gの二次勾配の平均値を利用して決定することができる。
　なお、ここでは具体的な実験データは示さないが、図３・図４の最適化アルゴリズムを用いることにより、従来の最適化アルゴリズム（例えば、図１・図２の最適化アルゴリズム）と比べて、行列{A, B}の収束速度、行列{A, B}の推定精度がともに向上することが確認された。特に、図３のB-P-R型の最適化アルゴリズムで効果が高いことが確認された。
＜第１実施形態＞
　Yを観測信号または観測データを表す非負の観測行列、Zを式(1-1)を満たす観測行列Yを構成する非負の行列とする。非負値行列分解最適化装置１００は、観測行列Yを入力として、Z=ABTを満たす行列Zの分解である非負の行列{A, B}を最適化する。
　以下、図５～図６を参照して非負行列分解最適化装置１００を説明する。図５は、非負行列分解最適化装置１００の構成を示すブロック図である。図６は、非負行列分解最適化装置１００の動作を示すフローチャートである。図５に示すように非負行列分解最適化装置１００は、最適化部１１０と、記録部１９０を含む。記録部１９０は、非負行列分解最適化装置１００の処理に必要な情報を適宜記録する構成部である。記録部１９０は、例えば、入力となる観測行列Y、最適化対象となる行列{A, B}を記録する。
　図６に従い非負行列分解最適化装置１００の動作について説明する。
　Ｓ１１０において、最適化部１１０は、観測行列Yを入力として、Z=ABTを満たす行列Zの分解である行列{A, B}を最適化する。具体的には、最適化部１１０は、観測行列Yを入力として、コスト関数Jの最小化問題を解くことにより行列{A, B}を最適化する。ここで、コスト関数J(A, B):RI×K×RJ×K→R∪{-∞, +∞}は、次式で定義される。
Figure JPOXMLDOC01-appb-M000066
（ただし、Gは損失項、HAは行列Aに対する正規化項、HBは行列Bに対する正規化項である。）
　以下、図７～図８を参照して最適化部１１０について説明する。図７は、最適化部１１０の構成を示すブロック図である。図８は、最適化部１１０の動作を示すフローチャートである。図７に示すように最適化部１１０は、初期化部１１１、行列Z更新部１１２と、第１双対変数更新部１１３と、行列A更新部１１４と、行列B更新部１１５と、第２双対変数更新部１１６と、カウンタ更新部１１７と、終了条件判定部１１８を含む。
　図８に従い最適化部１１０の動作について説明する。
　Ｓ１１１において、初期化部１１１は、カウンタtを初期化する。具体的には、t=1とする。また、初期化部１１１は、行列A、行列B、行列Z、双対変数~Xを初期化する。
　Ｓ１１２において、行列Z更新部１１２は、次式により、行列Zを更新する。
Figure JPOXMLDOC01-appb-M000067
（ただし、HZは行列Zに対する正規化項、BD^+は関数D+を用いて定義されるブレグマンダイバージェンスである。）
　なお、ブレグマンダイバージェンスBD^+の定義に用いる関数D+は、式
Figure JPOXMLDOC01-appb-M000068
（ただし、Cはcijを(i, j)要素とする行列、hij=∇ij 2G(cij old)（cij oldは現在の行列Cの(i, j) 要素、ε(>0)は所定の定数）により与えられる。また、関数D+は、hijを(i, j)要素とする行列Hに対してH=hAhB Tを満たすベクトルhA, hBを用いて計算してもよい。関数D+を計算する際、Zや~Xが用いられる。
　この関数D+は、行列A更新部１１４、行列B更新部１１５でも用いる。この場合、関数D+を計算する際、ABTや-~Vが用いられる。
　Ｓ１１３において、第１双対変数更新部１１３は、次式により、双対変数~Vを更新する。
Figure JPOXMLDOC01-appb-M000069
　Ｓ１１４において、行列A更新部１１４は、次式により、行列Aを更新する。
Figure JPOXMLDOC01-appb-M000070
　Ｓ１１５において、行列B更新部１１５は、次式により、行列Bを更新する。
Figure JPOXMLDOC01-appb-M000071
　Ｓ１１６において、第２双対変数更新部１１６は、次式により、双対変数~Xを更新する。
Figure JPOXMLDOC01-appb-M000072
　Ｓ１１７において、カウンタ更新部１１７は、カウンタtを1だけインクリメントする。具体的には、t←t+1とする。
　Ｓ１１８において、終了条件判定部１１８は、カウンタtが所定の更新回数T（Tは1以上の整数であり、例えば10万）に達した場合（つまり、t>Tとなり、終了条件が満たされた場合）は、そのときの行列{A, B}を出力して、処理を終了する。それ以外の場合、Ｓ１１２の処理に戻る。つまり、最適化部１１０は、Ｓ１１２～Ｓ１１８の処理を繰り返す。
　なお、Ｓ１１４とＳ１１５（つまり、行列Aの更新処理と行列Bの更新処理）については、Ｓ１１４、Ｓ１１５の順に実行してもよいし、Ｓ１１５、Ｓ１１４の順に実行してもよい。また、Ｓ１１４とＳ１１５を並列に実行するようにしてもよい。
（変形例）
　第２双対変数更新部１１６では、B-P-R型更新則により双対変数~Xを更新したが、B-D-R型更新則により更新するようにしてもよい。この場合、Ｓ１１６において、第２双対変数更新部１１６は、次式により、双対変数~Xを更新する。
Figure JPOXMLDOC01-appb-M000073
（ただし、ξは0<ξ<1を満たす定数）
　以上、まとめると、最適化部１１０は、ブレグマン単調作用素分解に基づいて行列{A, B}を最適化するものである。
　本実施形態の発明によれば、非負値行列を２つの非負値行列の積として表す非負値行列分解を高速かつ安定的に得ることができる。また、凸関数Gに合わせて関数D+を適応的に修正するため、変数空間の計量を適応的に修正することができ、従来のADMMに基づく最適化アルゴリズムで必要であった（例えば、ステップサイズなどの）パラメータの調整が不要となる。
＜補記＞
　本発明の装置は、例えば単一のハードウェアエンティティとして、キーボードなどが接続可能な入力部、液晶ディスプレイなどが接続可能な出力部、ハードウェアエンティティの外部に通信可能な通信装置（例えば通信ケーブル）が接続可能な通信部、ＣＰＵ（Central Processing Unit、キャッシュメモリやレジスタなどを備えていてもよい）、メモリであるＲＡＭやＲＯＭ、ハードディスクである外部記憶装置並びにこれらの入力部、出力部、通信部、ＣＰＵ、ＲＡＭ、ＲＯＭ、外部記憶装置の間のデータのやり取りが可能なように接続するバスを有している。また必要に応じて、ハードウェアエンティティに、ＣＤ－ＲＯＭなどの記録媒体を読み書きできる装置（ドライブ）などを設けることとしてもよい。このようなハードウェア資源を備えた物理的実体としては、汎用コンピュータなどがある。
　ハードウェアエンティティの外部記憶装置には、上述の機能を実現するために必要となるプログラムおよびこのプログラムの処理において必要となるデータなどが記憶されている（外部記憶装置に限らず、例えばプログラムを読み出し専用記憶装置であるＲＯＭに記憶させておくこととしてもよい）。また、これらのプログラムの処理によって得られるデータなどは、ＲＡＭや外部記憶装置などに適宜に記憶される。
　ハードウェアエンティティでは、外部記憶装置（あるいはＲＯＭなど）に記憶された各プログラムとこの各プログラムの処理に必要なデータが必要に応じてメモリに読み込まれて、適宜にＣＰＵで解釈実行・処理される。その結果、ＣＰＵが所定の機能（上記、…部、…手段などと表した各構成要件）を実現する。
　本発明は上述の実施形態に限定されるものではなく、本発明の趣旨を逸脱しない範囲で適宜変更が可能である。また、上記実施形態において説明した処理は、記載の順に従って時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されるとしてもよい。
　既述のように、上記実施形態において説明したハードウェアエンティティ（本発明の装置）における処理機能をコンピュータによって実現する場合、ハードウェアエンティティが有すべき機能の処理内容はプログラムによって記述される。そして、このプログラムをコンピュータで実行することにより、上記ハードウェアエンティティにおける処理機能がコンピュータ上で実現される。
　この処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体としては、例えば、磁気記録装置、光ディスク、光磁気記録媒体、半導体メモリ等どのようなものでもよい。具体的には、例えば、磁気記録装置として、ハードディスク装置、フレキシブルディスク、磁気テープ等を、光ディスクとして、ＤＶＤ（Digital Versatile Disc）、ＤＶＤ－ＲＡＭ（Random Access Memory）、ＣＤ－ＲＯＭ（Compact Disc Read Only Memory）、ＣＤ－Ｒ（Recordable）／ＲＷ（ReWritable）等を、光磁気記録媒体として、ＭＯ（Magneto-Optical disc）等を、半導体メモリとしてＥＥＰ－ＲＯＭ（Electronically Erasable and Programmable-Read Only Memory）等を用いることができる。
　また、このプログラムの流通は、例えば、そのプログラムを記録したＤＶＤ、ＣＤ－ＲＯＭ等の可搬型記録媒体を販売、譲渡、貸与等することによって行う。さらに、このプログラムをサーバコンピュータの記憶装置に格納しておき、ネットワークを介して、サーバコンピュータから他のコンピュータにそのプログラムを転送することにより、このプログラムを流通させる構成としてもよい。
　このようなプログラムを実行するコンピュータは、例えば、まず、可搬型記録媒体に記録されたプログラムもしくはサーバコンピュータから転送されたプログラムを、一旦、自己の記憶装置に格納する。そして、処理の実行時、このコンピュータは、自己の記憶装置に格納されたプログラムを読み取り、読み取ったプログラムに従った処理を実行する。また、このプログラムの別の実行形態として、コンピュータが可搬型記録媒体から直接プログラムを読み取り、そのプログラムに従った処理を実行することとしてもよく、さらに、このコンピュータにサーバコンピュータからプログラムが転送されるたびに、逐次、受け取ったプログラムに従った処理を実行することとしてもよい。また、サーバコンピュータから、このコンピュータへのプログラムの転送は行わず、その実行指示と結果取得のみによって処理機能を実現する、いわゆるＡＳＰ（Application Service Provider）型のサービスによって、上述の処理を実行する構成としてもよい。なお、本形態におけるプログラムには、電子計算機による処理の用に供する情報であってプログラムに準ずるもの（コンピュータに対する直接の指令ではないがコンピュータの処理を規定する性質を有するデータ等）を含むものとする。
　また、この形態では、コンピュータ上で所定のプログラムを実行させることにより、ハードウェアエンティティを構成することとしたが、これらの処理内容の少なくとも一部をハードウェア的に実現することとしてもよい。
　上述の本発明の実施形態の記載は、例証と記載の目的で提示されたものである。網羅的であるという意思はなく、開示された厳密な形式に発明を限定する意思もない。変形やバリエーションは上述の教示から可能である。実施形態は、本発明の原理の最も良い例証を提供するために、そして、この分野の当業者が、熟考された実際の使用に適するように本発明を色々な実施形態で、また、色々な変形を付加して利用できるようにするために、選ばれて表現されたものである。すべてのそのような変形やバリエーションは、公正に合法的に公平に与えられる幅にしたがって解釈された添付の請求項によって定められた本発明のスコープ内である。

Claims (7)

  1. 　Yを観測信号または観測データを表す非負の観測行列、Zを観測行列Yを構成する非負の行列とし、
    　観測行列Yを入力として、Z=ABTを満たす行列Zの分解である非負の行列{A, B}を最適化する最適化部を含む非負値行列分解最適化装置であって、
    　前記最適化部が最適化に用いるコスト関数Jは、式
    Figure JPOXMLDOC01-appb-M000001
     （ただし、Gは損失項、HAは行列Aに対する正規化項、HBは行列Bに対する正規化項）で定義されるものであり、
    　前記最適化部は、
    　次式により、行列Zを更新する行列Z更新部と、
    Figure JPOXMLDOC01-appb-M000002
     （ただし、HZは行列Zに対する正規化項、BD^+は関数D+を用いて定義されるブレグマンダイバージェンス）
    　次式により、双対変数~Vを更新する第１双対変数更新部と、
    Figure JPOXMLDOC01-appb-M000003
     　次式により、行列Aを更新する行列A更新部と、
    Figure JPOXMLDOC01-appb-M000004
     　次式により、行列Bを更新する行列B更新部と、
    Figure JPOXMLDOC01-appb-M000005
     　次式により、双対変数~Xを更新する第２双対変数更新部と、
    Figure JPOXMLDOC01-appb-M000006
     　を含む非負値行列分解最適化装置。
  2. 　Yを観測信号または観測データを表す非負の観測行列、Zを観測行列Yを構成する非負の行列とし、
    　観測行列Yを入力として、Z=ABTを満たす行列Zの分解である非負の行列{A, B}を最適化する最適化部を含む非負値行列分解最適化装置であって、
    　前記最適化部が最適化に用いるコスト関数Jは、式
    Figure JPOXMLDOC01-appb-M000007
     （ただし、Gは損失項、HAは行列Aに対する正規化項、HBは行列Bに対する正規化項）で定義されるものであり、
    　前記最適化部は、
    　次式により、行列Zを更新する行列Z更新部と、
    Figure JPOXMLDOC01-appb-M000008
     （ただし、HZは行列Zに対する正規化項、BD^+は関数D+を用いて定義されるブレグマンダイバージェンス）
    　次式により、双対変数~Vを更新する第１双対変数更新部と、
    Figure JPOXMLDOC01-appb-M000009
     　次式により、行列Aを更新する行列A更新部と、
    Figure JPOXMLDOC01-appb-M000010
     　次式により、行列Bを更新する行列B更新部と、
    Figure JPOXMLDOC01-appb-M000011
     　次式により、双対変数~Xを更新する第２双対変数更新部と、
    Figure JPOXMLDOC01-appb-M000012
     （ただし、ξは0<ξ<1を満たす定数）
    　を含む非負値行列分解最適化装置。
  3. 　請求項１または２に記載の非負値行列分解最適化装置であって、
    　前記行列Z更新部、前記行列A更新部、前記行列B更新部で用いるブレグマンダイバージェンスBD^+の定義に用いる関数D+は、式
    Figure JPOXMLDOC01-appb-M000013
     （ただし、Cはcijを(i, j)要素とする行列、hij=∇ij 2G(cij old)（cij oldは現在の行列Cの(i, j) 要素、ε(>0)は所定の定数）
    により与えられる
    　ことを特徴とする非負値行列分解最適化装置。
  4. 　請求項３に記載の非負値行列分解最適化装置であって、
    　関数D+は、hijを(i, j)要素とする行列Hに対してH=hAhB Tを満たすベクトルhA, hBを用いて計算される
    　ことを特徴とする非負値行列分解最適化装置。
  5. 　Yを観測信号または観測データを表す非負の観測行列、Zを観測行列Yを構成する非負の行列とし、
    　非負値行列分解最適化装置が、観測行列Yを入力として、Z=ABTを満たす行列Zの分解である非負の行列{A, B}を最適化する最適化ステップを実行する非負値行列分解最適化方法であって、
    　前記最適化ステップで最適化に用いるコスト関数Jは、式
    Figure JPOXMLDOC01-appb-M000014
     （ただし、Gは損失項、HAは行列Aに対する正規化項、HBは行列Bに対する正規化項）で定義されるものであり、
    　前記最適化ステップは、
    　次式により、行列Zを更新する行列Z更新ステップと、
    Figure JPOXMLDOC01-appb-M000015
     （ただし、HZは行列Zに対する正規化項、BD^+は関数D+を用いて定義されるブレグマンダイバージェンス）
    　次式により、双対変数~Vを更新する第１双対変数更新ステップと、
    Figure JPOXMLDOC01-appb-M000016
     　次式により、行列Aを更新する行列A更新ステップと、
    Figure JPOXMLDOC01-appb-M000017
     　次式により、行列Bを更新する行列B更新ステップと、
    Figure JPOXMLDOC01-appb-M000018
     　次式により、双対変数~Xを更新する第２双対変数更新ステップと、
    Figure JPOXMLDOC01-appb-M000019
     　を含む非負値行列分解最適化方法。
  6. 　Yを観測信号または観測データを表す非負の観測行列、Zを観測行列Yを構成する非負の行列とし、
    　非負値行列分解最適化装置が、観測行列Yを入力として、Z=ABTを満たす行列Zの分解である非負の行列{A, B}を最適化する最適化ステップを実行する非負値行列分解最適化方法であって、
    　前記最適化ステップで最適化に用いるコスト関数Jは、式
    Figure JPOXMLDOC01-appb-M000020
     （ただし、Gは損失項、HAは行列Aに対する正規化項、HBは行列Bに対する正規化項）で定義されるものであり、
    　前記最適化ステップは、
    　次式により、行列Zを更新する行列Z更新ステップと、
    Figure JPOXMLDOC01-appb-M000021
     （ただし、HZは行列Zに対する正規化項、BD^+は関数D+を用いて定義されるブレグマンダイバージェンス）
    　次式により、双対変数~Vを更新する第１双対変数更新ステップと、
    Figure JPOXMLDOC01-appb-M000022
     　次式により、行列Aを更新する行列A更新ステップと、
    Figure JPOXMLDOC01-appb-M000023
     　次式により、行列Bを更新する行列B更新ステップと、
    Figure JPOXMLDOC01-appb-M000024
     　次式により、双対変数~Xを更新する第２双対変数更新ステップと、
    Figure JPOXMLDOC01-appb-M000025
     （ただし、ξは0<ξ<1を満たす定数）
    　を含む非負値行列分解最適化方法。
  7. 　請求項１ないし４のいずれか１項に記載の非負値行列分解最適化装置としてコンピュータを機能させるためのプログラム。
PCT/JP2020/002206 2019-02-05 2020-01-23 非負値行列分解最適化装置、非負値行列分解最適化方法、プログラム Ceased WO2020162189A1 (ja)

Priority Applications (1)

Application Number Priority Date Filing Date Title
US17/428,263 US20220121728A1 (en) 2019-02-05 2020-01-23 Nonnegative matrix factorization optimization apparatus, nonnegative matrix factorization optimization method, and program

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
JP2019018425A JP7167746B2 (ja) 2019-02-05 2019-02-05 非負値行列分解最適化装置、非負値行列分解最適化方法、プログラム
JP2019-018425 2019-02-05

Publications (1)

Publication Number Publication Date
WO2020162189A1 true WO2020162189A1 (ja) 2020-08-13

Family

ID=71947810

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2020/002206 Ceased WO2020162189A1 (ja) 2019-02-05 2020-01-23 非負値行列分解最適化装置、非負値行列分解最適化方法、プログラム

Country Status (3)

Country Link
US (1) US20220121728A1 (ja)
JP (1) JP7167746B2 (ja)
WO (1) WO2020162189A1 (ja)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2022269680A1 (ja) * 2021-06-21 2022-12-29 日本電信電話株式会社 変数最適化システム
WO2024013865A1 (ja) * 2022-07-13 2024-01-18 日本電信電話株式会社 最適化問題求解装置、最適化問題求解方法、およびプログラム

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2023100221A1 (ja) * 2021-11-30 2023-06-08 日本電信電話株式会社 映像音声変換装置、方法およびプログラム

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2014203459A (ja) * 2013-04-03 2014-10-27 三菱電機株式会社 ２つ以上の光源の組によって照光されたシーンの２つ以上の入力画像の組から照明基底画像の組を復元する方法

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP5942420B2 (ja) * 2011-07-07 2016-06-29 ヤマハ株式会社 音響処理装置および音響処理方法
CN105989851B (zh) * 2015-02-15 2021-05-07 杜比实验室特许公司 音频源分离
JP7235129B2 (ja) * 2019-09-19 2023-03-08 日本電信電話株式会社 変数最適化装置、変数最適化方法、プログラム

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2014203459A (ja) * 2013-04-03 2014-10-27 三菱電機株式会社 ２つ以上の光源の組によって照光されたシーンの２つ以上の入力画像の組から照明基底画像の組を復元する方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
CHRETIEN, STEPHANE ET AL.: "ESTIMATING FEATURES WITH MISSING VALUES AND OUTLIERS: A BREGMAN-PROXIMAL POINT ALGORITHM FOR ROBUST NON-NEGATIVE MATRIX FACTORIZATION WITH APPLICATION TO GENE EXPRESSION ANALYSIS", ARXIV.ORG, 9 February 2015 (2015-02-09), XP021264558, Retrieved from the Internet <URL:https://arxiv.org/abs/1502.02523> *

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2022269680A1 (ja) * 2021-06-21 2022-12-29 日本電信電話株式会社 変数最適化システム
JPWO2022269680A1 (ja) * 2021-06-21 2022-12-29
JP7533788B2 (ja) 2021-06-21 2024-08-14 日本電信電話株式会社 変数最適化システム
WO2024013865A1 (ja) * 2022-07-13 2024-01-18 日本電信電話株式会社 最適化問題求解装置、最適化問題求解方法、およびプログラム
JPWO2024013865A1 (ja) * 2022-07-13 2024-01-18
JP7776788B2 (ja) 2022-07-13 2025-11-27 Ｎｔｔ株式会社 最適化問題求解装置、最適化問題求解方法、およびプログラム

Also Published As

Publication number Publication date
US20220121728A1 (en) 2022-04-21
JP2020126441A (ja) 2020-08-20
JP7167746B2 (ja) 2022-11-09

Similar Documents

Publication Publication Date Title
Ma et al. The Barron space and the flow-induced function spaces for neural network models
Adam et al. Image denoising using combined higher order non-convex total variation with overlapping group sparsity
Koblents et al. A population Monte Carlo scheme with transformed weights and its application to stochastic kinetic models
Landau et al. A polynomial time algorithm for the ground state of one-dimensional gapped local Hamiltonians
Leobacher et al. A numerical method for SDEs with discontinuous drift
WO2020162189A1 (ja) 非負値行列分解最適化装置、非負値行列分解最適化方法、プログラム
Prakash et al. Constrained nonlinear state estimation using ensemble Kalman filters
Castella et al. Stroboscopic averaging for the nonlinear Schrödinger equation
Li et al. Numerical solutions for weakly singular Volterra integral equations using Chebyshev and Legendre pseudo-spectral Galerkin methods
Hornik et al. On maximum likelihood estimation of the concentration parameter of von Mises–Fisher distributions
Ha et al. Resource analysis of quantum computing with noisy qubits for Shor’s factoring algorithms: J. Ha et al.
Bondarenko Direct and inverse problems for the matrix Sturm–Liouville operator with general self-adjoint boundary conditions
Bui et al. Gaussian process for predicting CPU utilization and its application to energy efficiency
Li et al. Finite element approximations of a class of nonlinear stochastic wave equations with multiplicative noise
Kuete Meli et al. A universal quantum algorithm for weighted maximum cut and Ising problems
JP7114497B2 (ja) 変数最適化装置、変数最適化方法、プログラム
Lam et al. Nonlinear shrinkage estimation of large integrated covariance matrices
Bader et al. Novel symplectic integrators for the Klein–Gordon equation with space-and time-dependent mass
Wu et al. A new non-convex sparse optimization method for image restoration
Cornebise et al. Adaptive sequential Monte Carlo by means of mixture of experts
Silva et al. An immune-inspired, information-theoretic framework for blind inversion of Wiener systems
JP7235129B2 (ja) 変数最適化装置、変数最適化方法、プログラム
Bakin et al. Parallel MARS algorithm based on B-splines
JP7156064B2 (ja) 潜在変数最適化装置、フィルタ係数最適化装置、潜在変数最適化方法、フィルタ係数最適化方法、プログラム
Zwaan et al. Multidirectional subspace expansion for one-parameter and multiparameter Tikhonov regularization

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: 20753159

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 20753159

Country of ref document: EP

Kind code of ref document: A1