WO2022113180A1 - 計算方法、計算装置およびプログラム - Google Patents

計算方法、計算装置およびプログラム Download PDF

Info

Publication number
WO2022113180A1
WO2022113180A1 PCT/JP2020/043722 JP2020043722W WO2022113180A1 WO 2022113180 A1 WO2022113180 A1 WO 2022113180A1 JP 2020043722 W JP2020043722 W JP 2020043722W WO 2022113180 A1 WO2022113180 A1 WO 2022113180A1
Authority
WO
WIPO (PCT)
Prior art keywords
square root
calculation result
calculation
computer
unit
Prior art date
Application number
PCT/JP2020/043722
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/JP2020/043722 priority Critical patent/WO2022113180A1/ja
Priority to JP2022564862A priority patent/JP7560765B2/ja
Publication of WO2022113180A1 publication Critical patent/WO2022113180A1/ja

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10Complex mathematical operations

Definitions

  • the present invention relates to a calculation method, a calculation device, and a program technology.
  • An embedded graph is defined by combining a graph, which is a set of a plurality of vertices and edges connecting the vertices, and the positions of the vertices.
  • a graph can be constructed by associating network nodes with vertices and network links with edges. At this time, if the physical positions of the network nodes and the node coordinates when visualizing the network correspond to the vertex positions, the embedded graph can be found.
  • borderline data such as borders, prefectural borders, coastlines, and lakeshore lines are generally held as embedded graphs.
  • the boundary line is expressed by plotting the vertices densely on the boundary line and connecting the sides.
  • embedded graphs are used to retain data related to shapes and regions. For example, data such as an object shape in a vector image, a depth boundary line in a depth map, a subject shape in a natural image, and a texture area shape are expressed as an embedded graph.
  • the boundaries are represented by edges.
  • the shape of the embedded graph is maintained as much as possible, and the number of vertices and sides is reduced as much as possible to simplify it (hereinafter, referred to as simplification of the embedded graph, or simply simplification).
  • simplification of the embedded graph or simply simplification.
  • the embedded graph for this network can be simplified, the network configuration can be presented in an easy-to-understand manner, the number of vertices and the number of sides can be reduced, data can be compressed, and the drawing load can be reduced.
  • the boundary line is accurately expressed by plotting the vertices densely on the boundary line, but the amount of data becomes enormous.
  • boundary line data By simplifying the boundary line data, it is possible to reduce the amount of data and transmit and store it, and it is also possible to generate visually easy-to-understand geographic information with simplified boundaries.
  • image engineering the larger the number of vertices and sides, the more detailed the data of the shape and region can be retained and accurately expressed.
  • the amount of data becomes enormous and the amount of code required for transmission and storage increases. It ends up.
  • the inverse square root a 1 / ⁇ A of a certain positive real number A is repeatedly calculated in the form of normal calculation in the reflected light and dispersed light simulation in the computer graphics (CG) field.
  • the high-speed reciprocal square root method (FastInverseSquareRoot; FastInvSqrt method) was invented, and by avoiding the calculation of the square root, the calculation of the reciprocal square root is speeded up several times.
  • This technique is implemented in game software using CG (see Non-Patent Document 2).
  • the modified FastInvSqrt method has been invented as a high-precision version of the FastInvSqrt method (see Non-Patent Document 3).
  • Japanese Unexamined Patent Publication No. 2016-221829 Japanese Unexamined Patent Publication No. 2017-211706
  • the graph simplification device using FastMultipleSVT cannot process at high speed, and graph drawing and graph processing cannot be executed at high speed.
  • the present invention aims to provide a technique for performing a calculation including a square root at a higher speed.
  • One aspect of the present invention is a first inverse square root step in which the computer calculates the inverse square root of the positive real number A by the fast inverse square root method, and the computer calculates the inverse square root of the positive real number B by the fast inverse square root method.
  • a second inverse square root step a subtraction step in which the computer subtracts B from A, and a computer multiplying the calculation result in the first inverse square root step by the calculation result in the second inverse square root step.
  • One aspect of the present invention is a program for causing a computer to execute the above calculation method.
  • One aspect of the present invention is a first inverse square root that calculates the inverse square root of a positive real number A by the fast inverse square root method, and a second inverse square root that calculates the inverse square root of the positive real B by the fast inverse square root method.
  • a subtraction part that subtracts B from A, a first multiplication part that multiplies the calculation result in the first inverse square root part, and the calculation result in the second inverse square root part, and the first multiplication part.
  • a second multiplication unit that multiplies the calculation result and the calculation result in the subtraction unit, and an addition unit that adds the calculation result in the first inverse square root portion and the calculation result in the second inverse square root portion.
  • a calculation device including a division unit that divides the calculation result in the second multiplication unit by the calculation result in the addition unit.
  • FIG. 1 is a diagram showing a configuration of an information processing device 1 including a calculation device 10 in the embodiment.
  • the information processing device 1 includes an input device 3, a calculation device 10, and an output device 5.
  • the input device 3 inputs input data such as numerical values to the calculation device 10.
  • the calculation device 10 performs a calculation using the input data, and outputs the calculation result to the output device 5.
  • the calculation device 10 includes a CPU (Central Processing Unit), a memory, an auxiliary storage device, etc. connected by a bus, and executes a calculation process by executing the calculation program. All or part of each function of the computing device 10 may be realized by using hardware such as ASIC (Application Specific Integrated Circuit), PLD (Programmable Logic Device), and FPGA (Field Programmable Gate Array).
  • the calculation program may be recorded on a computer-readable recording medium.
  • Computer-readable recording media include, for example, flexible disks, magneto-optical disks, ROMs, CD-ROMs, portable media such as semiconductor storage devices (for example, SSD: Solid State Drive), and storage of hard disks built in computer systems. It is a device.
  • the calculator may be transmitted over a telecommunication line.
  • the "high-speed reciprocal square root calculation unit" in each embodiment calculates the reciprocal square root by the high-speed reciprocal square root method.
  • FIG. 2 is a block diagram showing the configuration of the calculation device 10A for calculating the reciprocal square root and the square root difference in the embodiment.
  • the calculation device 10A includes a high-speed reciprocal square root calculation unit 11, 21, a subtraction unit 12, a multiplication unit 17, 117, an addition unit 13, and a division unit 14.
  • A is input to the high-speed reciprocal square root calculation unit 11 and the subtraction unit 12.
  • B is input to the high-speed reciprocal square root calculation unit 21 and the subtraction unit 12.
  • the subtraction unit 12 outputs AB to the multiplication unit 117.
  • the multiplication unit 17 outputs ab to the multiplication unit 117.
  • the multiplication unit 117 outputs (AB) ab to the division unit 14.
  • the addition unit 13 outputs a + b to the division unit 14.
  • the division unit 14 outputs (AB) ab / (a + b) as c to the output device 5.
  • the arithmetic unit 10A calculates the reciprocal square root, but does not calculate ⁇ A and ⁇ B. Therefore, c is output without calculating ⁇ A and ⁇ B. Further, by rationalizing the numerator, it is possible to prevent the digit loss of ⁇ A- ⁇ B as compared with the case of calculating ⁇ A and ⁇ B.
  • FIG. 3 is a block diagram showing a calculation device 20 that calculates ⁇ A and ⁇ B and outputs a, b, and c as in the conventional case.
  • the calculation device 20 includes square root calculation units 15 and 25, a subtraction unit 22, and a division unit 24 and 34.
  • the square root calculation unit 15 When a positive real number A is input in the arithmetic unit 20, the square root calculation unit 15 outputs ⁇ A to the subtraction unit 22 and the division unit 24. When a positive real number B is input, the square root calculation unit 25 outputs ⁇ B to the subtraction unit 22 and the division unit 34.
  • FIG. 4 is a diagram showing the time required to calculate each calculation 1,000,000 times in C ++ language, with the addition being 1 flops.
  • sqrt (x) shows the square root
  • 1 / sqrt (x) shows the amount of calculation when the reciprocal square root is calculated by the standard library function of the C ++ language.
  • fastInvSqrt (x) shows the amount of calculation by the Fast Inverse Square Root method once in Newton's method.
  • fastInvSqrt2 (x) shows the amount of calculation by the Fast Inverse Square Root method twice in Newton's method.
  • fastInvSqrt3 (x) shows the amount of calculation by the Fast Inverse Square Root method three times in Newton's method.
  • mFastInvSqrt (x) shows the amount of calculation by the modified Fast Inverse Square Root method.
  • the total amount of calculation is 17 flops.
  • the calculation device 10A according to the present embodiment is used for calculation, the total amount of calculation is 10 flops. Therefore, it can be estimated that the arithmetic unit 10A according to the present embodiment reduces the floating-point arithmetic of 7 flops by about 41.2% as compared with the conventional arithmetic unit 20.
  • FIG. 5 is a block diagram showing a configuration of the calculation device 10B for calculating the reciprocal square root of the sum and the difference and the square root difference of the sum and the difference in the embodiment.
  • the arithmetic unit 10B includes a high-speed reciprocal square root calculation unit 31, 41, a subtraction unit 32, a multiplication unit 27, 37, 47, an addition unit 23, 33, and a division unit 44.
  • X When a positive real number X is input in the arithmetic unit 10B, X is input to the addition unit 23 and the subtraction unit 32.
  • e When a positive real number e is input, e is input to the addition unit 23, the subtraction unit 32, and the multiplication unit 27.
  • the addition unit 23 outputs X + e to the high-speed reciprocal square root calculation unit 31.
  • the subtraction unit 32 outputs X-e to the high-speed reciprocal square root calculation unit 41.
  • the multiplication unit 27 outputs 2e to the multiplication unit 47.
  • the multiplication unit 47 outputs 2eab to the division unit 44.
  • the division unit 44 outputs 2eab / (a + b) as c to the output device 5.
  • 2eab / (a + b) ⁇ (X + e) - ⁇ (X-e) is shown by letting A be X + e and B be X-e in the first embodiment.
  • the arithmetic unit 10B calculates the reciprocal square root, but does not calculate ⁇ (X + e) and ⁇ (X-e). Therefore, c is output without calculating ⁇ (X + e) and ⁇ (X-e).
  • FIG. 6 is a block diagram showing a calculation device 40 that calculates ⁇ (X + e) and ⁇ (X-e) and outputs a, b, and c as in the conventional case.
  • the calculation device 40 includes a square root calculation unit 35, 45, a subtraction unit 42, 52, an addition unit 43, and a division unit 54, 64.
  • the addition unit 43 When the positive real numbers X and e are input in the calculation device 40, the addition unit 43 outputs X + e to the square root calculation unit 35.
  • the subtraction unit 42 outputs X-e to the square root calculation unit 45.
  • the square root calculation unit 35 outputs ⁇ (X + e) to the division unit 54 and the subtraction unit 52.
  • the square root calculation unit 45 outputs ⁇ (X—e) to the division unit 64 and the subtraction unit 52.
  • the total amount of calculation is 19 flops.
  • the calculation device 10B according to the present embodiment is used for calculation, the total amount of calculation is 12 flops. Therefore, it can be estimated that the arithmetic unit 10B according to the present embodiment reduces the floating-point arithmetic of 7 flops by about 36.8% as compared with the conventional arithmetic unit 40.
  • FIG. 7 is a block diagram showing the configuration of the calculation device 10C for calculating the second singular value of the M ⁇ 2 matrix, the reciprocal of the sum of the singular values, and the reciprocal of the difference of the singular values in the embodiment.
  • It is a calculation device that outputs the reciprocal 1 / ( ⁇ 1 ⁇ ⁇ 2 ) of the difference between the two to the output device 5.
  • the arithmetic unit 10C includes the arithmetic unit 10B described in the second embodiment. Further, the arithmetic unit 10C includes a decomposition unit 110, an inner product unit 16, 26, 36, a subtraction unit 62, a multiplication unit 57, 67, 77, 87, an addition unit 53, and a square root calculation unit 55.
  • Y is input to the decomposition unit 110.
  • the decomposition unit 110 decomposes the M ⁇ 2 matrix Y into the column vectors y 1 and y 2 .
  • the decomposition unit 110 outputs y 1 to the inner product units 16 and 36.
  • the decomposition unit 110 outputs y 2 to the inner product units 26 and 36.
  • the inner product unit 16 calculates the inner product a of y 1 and y 1 , and outputs a to the addition unit 53, the multiplication unit 67, and the output device 5.
  • the inner product unit 26 calculates the inner product c of y 2 and y 2 , and outputs c to the addition unit 53, the multiplication unit 67, and the output device 5.
  • the inner product unit 36 calculates the inner product b of y 1 and y 2 , and outputs b to the multiplication unit 77 and the output device 5.
  • the multiplication unit 67 outputs ac to the subtraction unit 62.
  • the multiplying unit 77 outputs b 2 to the subtracting unit 62.
  • the multiplying unit 57 outputs 2e to the arithmetic unit 10B.
  • the calculation device 10B outputs 1 / ⁇ (f + 2e) to the output device 5 as 1 / ( ⁇ 1 + ⁇ 2 ).
  • the calculation device 10B outputs 1 / ⁇ (f-2e) to the output device 5 as 1 / ( ⁇ 1 ⁇ ⁇ 2 ).
  • the arithmetic unit 10B outputs ⁇ (f + 2e) ⁇ (f-2e) as 2 ⁇ 2 to the multiplication unit 87.
  • the multiplying unit 87 outputs ⁇ 2 to the output device 5.
  • 1 / ⁇ (f ⁇ 2e) is 1 / ( ⁇ 1 ⁇ ⁇ 2 ) as follows.
  • the computer 10C has a second singular value ⁇ 2 , a reciprocal of the sum of the singular values 1 / ( ⁇ 1 + ⁇ 2 ), and a reciprocal of the difference between the singular values 1 / ( ⁇ 1 ⁇ ⁇ ). Not only 2 ) but also a, b, c, d, e, and f are output to the output device 5.
  • the calculation amount of the calculation device 10C will be described. Floating-point arithmetic is not performed (0 flop) because the processing of the decomposition unit 110 only decomposes the matrix into column vectors. Further, the amount of calculation by the inner product portions 16, 26, 36 is the amount of calculation of 2M-1 flops with respect to the input of the M-th order column vector. The square root calculation by sqrt is 4 flops as described above. Therefore, the total amount of calculation by the arithmetic unit 10C is 6M + 19 flops. As described above, since the arithmetic unit 10B has the effect of reducing 7 drops, the arithmetic unit 10C also has the effect of reducing 7 drops by the arithmetic unit 10B.
  • FIG. 8 is a block diagram showing the configuration of the calculation device 10D for calculating the second singular value of the 2 ⁇ 2 matrix, the reciprocal of the sum of the singular values, and the reciprocal of the difference of the singular values in the embodiment.
  • a 2 ⁇ 2 matrix Y [y 1 , y 2 ] is input from the input device 3, the second singular value ⁇ 2 , the reciprocal of the sum of the singular values 1 / ( ⁇ 1 + ⁇ 2 ), and the singular value.
  • It is a calculation device that outputs the reciprocal 1 / ( ⁇ 1 ⁇ ⁇ 2 ) of the difference between the two to the output device 5.
  • the arithmetic unit 10D includes the arithmetic unit 10B described in the second embodiment. Further, the arithmetic unit 10D includes a decomposition unit 120, an inner product unit 46, 56, a determinant calculation unit 19, an absolute value calculation unit 18, a multiplication unit 97, 107, and an addition unit 63.
  • Y is input to the decomposition unit 120.
  • the decomposition unit 120 decomposes the 2 ⁇ 2 matrix Y into the column vectors y 1 and y 2 .
  • the decomposition unit 120 outputs y 1 to the inner product unit 46.
  • the decomposition unit 120 outputs y 2 to the inner product unit 56.
  • the inner product unit 46 calculates the inner product a of y 1 and y 1 , and outputs a to the addition unit 63.
  • the inner product unit 56 calculates the inner product c of y 2 and y 2 , and outputs c to the addition unit 63.
  • the determinant calculation unit 19 calculates the determinant d of Y and outputs d to the absolute value calculation unit 18 and the output device 5.
  • the absolute value calculation unit 18 outputs the absolute value e of d to the multiplication unit 97.
  • the multiplying unit 97 outputs 2e to the arithmetic unit 10B.
  • the calculation device 10B outputs 1 / ⁇ (f + 2e) to the output device 5 as 1 / ( ⁇ 1 + ⁇ 2 ).
  • the calculation device 10B outputs 1 / ⁇ (f-2e) to the output device 5 as 1 / ( ⁇ 1 ⁇ ⁇ 2 ).
  • the arithmetic unit 10B outputs ⁇ (f + 2e) ⁇ (f-2e) as 2 ⁇ 2 to the multiplication unit 107.
  • the multiplying unit 107 outputs ⁇ 2 to the output device 5.
  • 1 / ⁇ (f ⁇ 2e) is 1 / ( ⁇ 1 ⁇ ⁇ 2 ) as shown in the third embodiment.
  • the computer 10D has a second singular value ⁇ 2 , a reciprocal of the sum of the singular values 1 / ( ⁇ 1 + ⁇ 2 ), and a reciprocal of the difference between the singular values 1 / ( ⁇ 1 ⁇ ⁇ ). Not only 2 ) but also d and f are output to the output device 5.
  • the calculation amount of the calculation device 10D will be described.
  • the amount of calculation of the determinant calculation unit 19 is 3 flops. Since the absolute value calculation unit 18 only evaluates the sign and inverts it when it is negative, the floating point operation is not performed (0 flop).
  • the total amount of calculation by the arithmetic unit 10D is 24 flops. As described above, since the arithmetic unit 10B has the effect of reducing 7 drops, the arithmetic unit 10D also has the effect of reducing 7 drops by the arithmetic unit 10B. [Fifth Embodiment] Hereinafter, the arithmetic unit 10E, which is the arithmetic unit 10 in the fifth embodiment, will be described. FIG.
  • FIG. 9 is a block diagram showing the configuration of the calculation device 10E for calculating the SVT of the M ⁇ 2 matrix in the embodiment.
  • SVT singular value thresholding
  • I 2 in (1) is a 2 ⁇ 2 identity matrix. Further, ⁇ , ⁇ , e, a, b, c, and d will be described later.
  • the calculation device 10E includes the calculation device 10C described in the third embodiment, the coefficient calculation device 200, and the M ⁇ 2 matrix conversion device 300.
  • the computing device 10C has a second singular value ⁇ 2 , a reciprocal of the sum of the singular values 1 / ( ⁇ 1 + ⁇ 2 ), and a reciprocal of the difference of the singular values 1 / ( ⁇ 1 ⁇ ⁇ ). 2 ), a, b, c, d, e, f are output.
  • d, f, ⁇ 2 , g, and h are output to the coefficient calculation device 200.
  • a, b, c, and e are output to the M ⁇ 2 matrix conversion device 300.
  • the coefficient calculation device 200 In the coefficient calculation device 200, d, f, ⁇ 2 , g, h, and ⁇ are input. The coefficient calculation device 200 outputs ⁇ and ⁇ .
  • the calculation method of ⁇ and ⁇ will be described. First, let R be a sufficiently large real number. Specifically, the size of R is about one tenth of the maximum number of single-precision floating-point types. Then, the coefficient calculation device 200 calculates ⁇ and ⁇ by performing the following four cases (a) to (d), and outputs them.
  • the reciprocal square root by the high-speed reciprocal square root method may be expressed as isqrt ( ⁇ ).
  • the reciprocal square root of a positive real number A calculated by the high-speed reciprocal square root method may be expressed as isqrt (A).
  • + in the lower right subscript in ( ⁇ ) + indicates a ramp function.
  • the ramp function outputs 0 if the input is a negative real number, and outputs the input real number as it is if the input is a non-negative real number.
  • the function min outputs the smallest value among the input numerical values.
  • the calculation amount of the calculation device 10E will be described.
  • Ramp function and function min are mainly calculated by comparing real values, and floating point arithmetic is not performed (0 flop).
  • the amount of calculation is 12M + 38flops.
  • the total amount of calculation by the arithmetic unit 10E is 12M + 29flops. From the above, it can be estimated that the floating-point arithmetic of 9 flops will be reduced.
  • the arithmetic unit 10F which is the arithmetic unit 10 in the sixth embodiment, will be described.
  • FIG. 10 is a block diagram showing a configuration of a calculation device 10F for calculating an SVT of a 2 ⁇ 2 matrix in an embodiment.
  • y ij is a component of Y in rows and columns of i.
  • the function sign ( ⁇ ) in (2) is a sign function, and outputs -1 when the input is negative, outputs 0 when the input is 0, and outputs +1 when the input is positive. .. Further, ⁇ , ⁇ , and d will be described later.
  • the calculation device 10F includes the calculation device 10D described in the fourth embodiment, the coefficient calculation device 201, and the 2 ⁇ 2 matrix conversion device 301.
  • d is output to the 2 ⁇ 2 matrix conversion device 300.
  • the coefficient calculation device 201 In the coefficient calculation device 201, d, f, ⁇ 2 , g, h, and ⁇ are input. The coefficient calculation device 201 outputs ⁇ and ⁇ . The calculation method of ⁇ and ⁇ will be described. First, let R be a sufficiently large real number. Specifically, the size of R is about one tenth of the maximum number of single-precision floating-point types. Then, the coefficient calculation device 200 calculates ⁇ and ⁇ by performing the following four cases (a) to (d), and outputs them.
  • the coefficient calculation device 200 outputs ⁇ and ⁇ to the 2 ⁇ 2 matrix conversion device 301 according to the above cases. As a result, since the 2 ⁇ 2 matrix conversion device 301 has all the parameters included in the above (2), (2) is output to the output device 5 using them.
  • the amount of calculation of the arithmetic unit 10F will be explained.
  • the amount of calculation is 41 flops.
  • the total amount of calculation by the arithmetic unit 10F is 32 flops. From the above, it can be estimated that the floating-point arithmetic of 9 flops will be reduced.
  • FIG. 11 is a diagram showing an outline of the graph simplification process.
  • G (V, E, P) in FIG. 11 indicates a graph embedded in the M-dimensional space.
  • V and E are sets of vertices and edges, respectively, and P ⁇ R
  • ⁇ M is the vertex coordinates (R here is a set of all real numbers).
  • ) T represents the coordinates of each vertex of the graph.
  • V ⁇ ⁇ v ⁇ V
  • degv 2 ⁇ be the set of vertices having a degree of 2.
  • FIG. 12 is a diagram showing a configuration example of the graph simplification device 500.
  • the graph simplification device 500 inputs the threshold value ⁇ and the graph G, and outputs the simplified graph G ′′.
  • the graph simplification device 500 includes a local linear alignment unit 510 and an unnecessary vertex removal unit 520.
  • the local linear alignment unit 510 outputs G'in which the graph G is locally linearly aligned using the threshold value ⁇ to the unnecessary vertex removal unit 520.
  • the unnecessary vertex removing unit 520 outputs G ′′ from which unnecessary vertices have been removed from the input vertices of the graph G ′.
  • FIG. 13 is a diagram showing a configuration example of the local linear alignment unit 510 in FIG. 12.
  • the local linear alignment unit 510 outputs G'in which the graph G is locally linearly aligned using the threshold value ⁇ to the unnecessary vertex removal unit 520.
  • the local linear alignment unit 510 includes a convex optimization problem formulating unit 511, a convex optimization problem solving unit 512, and a coordinate information replacement unit 513.
  • FIG. 14 is a diagram for explaining an operation of arranging vectors of adjacent sides.
  • the matrix LV shown in FIG. 14 is a matrix for arranging vectors from v to adjacent vertices. For the elements of the matrix LV shown in FIG. 14, all the “...” parts are 0. Further, when the coordinates of v are the k-th row of P , "-1" in LV shown in FIG. 14 is the k-th column. LVP gives a vector starting from v and ending at the coordinates of the k-1th row and a vector starting from v and ending at the coordinates of the k + 1th row. The operation that makes this LV act is a linear map.
  • the scale of faithful reproduction is the L1 norm
  • the regularization of the number of bends of the sides is the nuclear norm function
  • the convex optimization problem formulating unit 511 formulates the optimization problem as described in (3) below.
  • the karyotype norm function is a function that calculates the sum of the singular values of the input matrix.
  • FIG. 15 shows a specific procedure executed by the convex optimization problem solving unit 512.
  • FIG. 15 is a diagram showing an algorithm showing a solution of a local linear alignment problem.
  • f in prox ⁇ f of the 4th line of FIG. 15 is the following (4).
  • the above (5) is an SVT with the threshold value ⁇ of the matrix of the following (6).
  • the processing on the 8th line occupies about 53.8% of the calculation time of the algorithm shown in FIG. Therefore, by calculating the SVT using the above-mentioned calculation devices 10E and 10F, the graph simplification process can be executed at a higher speed than in the conventional case.
  • the coordinate information replacement unit 513 replaces the coordinate information with the vertex coordinates locally linearly aligned by the convex optimization problem solving unit 512, and outputs the locally linearly aligned G'.
  • This embodiment can be applied not only to the simplification of the graph but also to all the processes for performing SVT. For example, image false color removal falls into the problem of regularizing a large number of small matrices, as well as graph simplification.
  • the algorithm can be easily parallelized to data, and a large number of matrices can be processed simultaneously. If the algorithm can be parallelized to data, the speed can be increased by implementing a data parallel architecture such as Single Instrument Multiple Data (SIMD) adopted by many CPUs installed in personal computers.
  • SIMD Single Instrument Multiple Data
  • a storage device such as a memory in which various calculations such as an addition unit and an inner product unit temporarily store the calculation results is provided, and the calculation results are temporarily stored in the storage device. You may.
  • the present invention is applicable to a calculation device that calculates a square root.
  • Coefficient calculation device 300, 301 ... Matrix conversion device, 500 ... Graph simplification device, 510 ... Local linear alignment unit, 511 ... Convex optimization problem formula unit, 512 ... Convex optimization problem solution unit, 513 ... Coordinate information replacement part 520 ... Unnecessary vertex removal part

Landscapes

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

Abstract

本発明の一態様は、コンピュータが、正の実数Aの逆数平方根を高速逆数平方根法により計算する第1逆数平方根ステップと、コンピュータが、正の実数Bの逆数平方根を高速逆数平方根法により計算する第2逆数平方根ステップと、コンピュータが、AからBを減算する減算ステップと、コンピュータが、第1逆数平方根ステップでの計算結果と、第2逆数平方根ステップでの計算結果とを乗算する第1乗算ステップと、コンピュータが、第1乗算ステップでの計算結果と、減算ステップでの計算結果とを乗算する第2乗算ステップと、コンピュータが、第1逆数平方根ステップでの計算結果と、第2逆数平方根ステップでの計算結果とを加算する加算ステップと、コンピュータが、第2乗算ステップでの計算結果を、加算ステップでの計算結果で除算する除算ステップと、を備えた計算方法である。

Description

計算方法、計算装置およびプログラム
 本発明は、計算方法、計算装置およびプログラムの技術に関する。
 ネットワークや地理情報学、画像工学などの諸分野において、埋め込みグラフというデータ構造が存在する。埋め込みグラフとは、複数の頂点とその頂点同士を結ぶ辺の集合であるグラフと、その頂点位置を合わせて定義したものである。例えばネットワークでは、ネットワークノードを頂点、ネットワークリンクを辺に対応させてグラフを構成できる。このとき、ネットワークノードの物理的な位置や、ネットワークを可視化する際のノード座標を頂点位置に対応させると、埋め込みグラフを見出すことができる。
 また、地理情報学であれば、国境や県境、海岸線、湖岸線等の境界線データは一般に、埋め込みグラフとしてデータが保持されている。境界線上に密に頂点をプロットして辺を結ぶことで、境界線を表現している。また、画像工学分野であれば、形状や領域に関するデータの保持に埋め込みグラフが使用される。例えばベクター画像におけるオブジェクト形状、深度マップにおける深度境界線、自然画像の被写体形状やテクスチャ領域形状等のデータは、埋め込みグラフとして表現される。形状や領域の境界上に頂点をプロットすることで、境界線を辺により表現している。
 以上で例示した埋め込みグラフにおいては、埋め込みグラフの形状をできるだけ保ちつつ、頂点数や辺数を可能な限り削減して単純化すること(以下、埋め込みグラフの単純化、あるいは単に単純化と呼ぶ)に大きなメリットがある。
 ネットワークにおいては、頂点数、辺数が多いほど詳細なデータを表現できる一方で、一見してネットワーク構成の概要を掴むことが困難になり、加えてデータ量が膨大になるため描画に大きな負荷がかかる。このネットワークについての埋め込みグラフを単純化できれば、ネットワーク構成を視覚的に分かりやすく提示でき、頂点数、辺数を減らしてデータ圧縮し、描画の負荷を低減することが可能になる。地理情報学においては、境界線上に密に頂点をプロットすることで正確に境界線を表現するが、データ量は膨大になってしまう。
 境界線データを単純化することで、データ量を削減して伝送や蓄積が可能な他、境界を簡略した視覚的に分かりやすい地理情報を生成することができる。また画像工学においては、頂点数、辺数が多いほど形状や領域のデータを詳細に保持して正確に表現可能である一方で、データ量が膨大になり、伝送や蓄積に要する符号量が増えてしまう。
 この形状や領域を単純化できれば、表現の正確性を可能な限り保ちながら、頂点数、辺数を削減して符号量を削減できる。
 以上の単純化を達成するための発明として、埋め込みグラフ単純化法が提案されている(特許文献2参照)。またこの高速処理法(2次元版)(特許文献3参照)と、高速処理法(多次元版)(特許文献4参照)とが提案されている。
 上記の高速処理法では特異値閾値処理(Singular Value Thresholding;SVT)(特許文献4参照)を多数並列に高速計算するFast Multiple SVT(非特許文献1参照)が採用されている。図16にアルゴリズムを掲載する。図16Aと図16Bに示されるアルゴリズムのうちの3行目g-1とh-1と4行目σにおいてはそれぞれ逆数平方根と平方根の差が計算されている。この逆数平方根は図16Aと図16Bに示されるアルゴリズムの計算速度のボトルネックとなっている。
 また平方根の差の計算は桁落ちが生じやすく、誤差が発生しやすい。さらに、図16Aと図16Bに示されるアルゴリズムのうちの4行目の平方根の差計算においてのみ、平方根の計算が必要であり、その他の個所では逆数平方根しか用いられていない。
 さて、ある正の実数Aの逆数平方根a=1/√Aは、コンピュータグラフィクス(CG)分野における反射光および分散光シミュレーションにおいて、法線計算という形で多量に繰り返して計算される。
 この計算量を削減するために、高速逆数平方根法(Fast Inverse Square Root;FastInvSqrt法)が発明され、平方根の計算を回避することで、逆数平方根の計算は数倍程度高速化されている。この技術はCGを用いたゲームソフトウェアに実装されている(非特許文献2参照)。またFastInvSqrt法の高精度版として修正FastInvSqrt法が発明されている(非特許文献3参照)。
 このFastInvSqrt法または修正FastInvSqrt法を、上述のFast Multiple SVTにおける逆数平方根の計算に採用することで計算速度の向上が見込まれる。
特開2016-221829号公報 特開2017-211706号公報 特開2018-082249号公報 特開2019-046196号公報
佐々木崇元,北原正樹,清水淳,"低ランク最適化のための高速特異値閾値処理の数理," 第16回情報科学技術フォーラム,2017 M. Robertson, "A Brief History of InvSqrt," Bachelor Thesis, University of New Brunswick, 2012. C. J. Walczyk, ''A Modification of the Fast Inverse Square, '' MDPI Computation, vol. 7, no. 3, 2019
 しかしながら、Fast Multiple SVTでは平方根の差計算も必要であるため、結局のところ平方根が必要であり、さらに逆数演算を取る必要が発生する。このため、上記の速度向上は達成されない。
 このように、Fast Multiple SVTを用いるグラフ単純化装置は高速に処理できず、グラフ描画やグラフ処理を高速に実行できない。
 上記事情に鑑み、本発明は、平方根を含む計算をより速く行う技術の提供を目的としている。
 本発明の一態様は、コンピュータが、正の実数Aの逆数平方根を高速逆数平方根法により計算する第1逆数平方根ステップと、コンピュータが、正の実数Bの逆数平方根を高速逆数平方根法により計算する第2逆数平方根ステップと、コンピュータが、AからBを減算する減算ステップと、コンピュータが、前記第1逆数平方根ステップでの計算結果と、前記第2逆数平方根ステップでの計算結果とを乗算する第1乗算ステップと、コンピュータが、前記第1乗算ステップでの計算結果と、前記減算ステップでの計算結果とを乗算する第2乗算ステップと、コンピュータが、前記第1逆数平方根ステップでの計算結果と、前記第2逆数平方根ステップでの計算結果とを加算する加算ステップと、コンピュータが、前記第2乗算ステップでの計算結果を、前記加算ステップでの計算結果で除算する除算ステップと、を備えた計算方法である。
 本発明の一態様は、上記の計算方法をコンピュータに実行させるためのプログラムである。
 本発明の一態様は、正の実数Aの逆数平方根を高速逆数平方根法により計算する第1逆数平方根部と、正の実数Bの逆数平方根を高速逆数平方根法により計算する第2逆数平方根部と、AからBを減算する減算部と、前記第1逆数平方根部での計算結果と、前記第2逆数平方根部での計算結果とを乗算する第1乗算部と、前記第1乗算部での計算結果と、前記減算部での計算結果とを乗算する第2乗算部と、前記第1逆数平方根部での計算結果と、前記第2逆数平方根部での計算結果とを加算する加算部と、前記第2乗算部での計算結果を、前記加算部での計算結果で除算する除算部と、を備えた計算装置である。
 本発明により、平方根を含む計算をより速く行うことが可能となる。
計算装置を含む情報処理装置の構成を示す図である。 計算装置10Aの構成を示すブロック図である。 計算装置20を示すブロック図である。 計算量を示す図である。 計算装置10Bの構成を示すブロック図である。 計算装置40を示すブロック図である。 計算装置10Cの構成を示すブロック図である。 計算装置10Dの構成を示すブロック図である。 計算装置10Eの構成を示すブロック図である。 計算装置10Fの構成を示すブロック図である。 グラフの単純化処理の処理概要を示す図である。 グラフ単純化装置の構成例を示す図である。 局所線形整列部の構成例を示す図である。 隣接辺のベクトルを並べる操作を説明するための図である。 局所線形整列問題の解法を示すアルゴリズムを示す図である。 局所線形整列問題の解法を示すアルゴリズムを示す図である。 特異値閾値処理のアルゴリズムを示す図である。 特異値閾値処理のアルゴリズムを示す図である。
 本発明の実施形態について、図面を参照して詳細に説明する。
 図1は、実施形態における計算装置10を含む情報処理装置1の構成を示す図である。情報処理装置1は、入力装置3、計算装置10、および出力装置5を備える。入力装置3は、計算装置10に数値等の入力データを入力する。計算装置10は、入力データを用いて計算を行い、計算結果を出力装置5に出力する。
 計算装置10は、バスで接続されたCPU(Central Processing Unit)やメモリや補助記憶装置などを備え、計算プログラムに実行することによって計算処理を実行する。計算装置10の各機能の全て又は一部は、ASIC(Application Specific Integrated Circuit)やPLD(Programmable Logic Device)やFPGA(Field Programmable Gate Array)等のハードウェアを用いて実現されてもよい。計算プログラムは、コンピュータ読み取り可能な記録媒体に記録されてもよい。コンピュータ読み取り可能な記録媒体とは、例えばフレキシブルディスク、光磁気ディスク、ROM、CD-ROM、半導体記憶装置(例えばSSD:Solid State Drive)等の可搬媒体、コンピュータシステムに内蔵されるハードディスク等の記憶装置である。計算プログラムは、電気通信回線を介して送信されてもよい。
 以下、計算装置10の第1実施形態~第6実施形態について説明する。その後、計算装置10の適用例としてグラフの単純化処理について説明する。なお、各実施形態における「高速逆数平方根計算部」は、高速逆数平方根法により逆数平方根を計算する。
 [第1実施形態]
 以下、第1実施形態における計算装置10である計算装置10Aについて説明する。図2は、実施形態における逆数平方根、および平方根差を計算する計算装置10Aの構成を示すブロック図である。計算装置10Aは、入力装置3から正の実数A、Bが入力され、a=1/√A、b=1/√B、c=√A-√Bを出力装置5に出力する計算装置である。計算装置10Aは、高速逆数平方根計算部11、21、減算部12、乗算部17、117、加算部13、および除算部14を備える。
 計算装置10Aにおいて、正の実数Aが入力されると、Aは、高速逆数平方根計算部11および減算部12に入力される。正の実数Bが入力されると、Bは、高速逆数平方根計算部21および減算部12に入力される。高速逆数平方根計算部11は、減算部12、加算部13、および出力装置5にa=1/√Aを出力する。高速逆数平方根計算部21は、減算部12、加算部13、および出力装置5にb=1/√Bを出力する。減算部12は、A-Bを乗算部117に出力する。
 乗算部17は、乗算部117にabを出力する。乗算部117は、除算部14に(A-B)abを出力する。加算部13は、a+bを除算部14に出力する。除算部14は、cとして(A-B)ab/(a+b)を出力装置5に出力する。
 ここで、(A-B)ab/(a+b)=√A-√Bであることは、以下の通りである。
√A-√B=(√A-√B)(√A+√B)/(√A+√B)
     =(A-B)/(√A+√B)
ここで、
√A+√B=(1/√A+1/√B)(√A√B)=(a+b)/(ab)
よって
(A-B)/(√A+√B)=(A-B)ab/(a+b)
 以上より、計算装置10Aは、a=1/√A、b=1/√B、およびc=√A-√Bを出力装置5に出力する。計算装置10Aは、逆数平方根を計算するが、√A、√Bは計算していない。よって、cを√A、√Bを計算することなく出力する。また、分子の有理化を行うことで、√A、√Bを計算する場合と比較して、√A-√Bの桁落ちを防止することができる。
 図3は、参考例として、従来通り、√A、√Bを計算してa,b,cを出力する計算装置20を示すブロック図である。計算装置20は、平方根計算部15、25、減算部22、および除算部24、34を備える。
 計算装置20において、正の実数Aが入力されると、平方根計算部15は、減算部22、および除算部24に√Aを出力する。正の実数Bが入力されると、平方根計算部25は、減算部22、および除算部34に√Bを出力する。
 除算部24は、a=1/√Aを出力する。除算部34は、b=1/√Bを出力する。減算部22は、c=√A-√Bを出力する。
 ここで計算量を比較するために、図4を用いて各計算の計算量について説明する。図4は、C++言語で、各計算を1,000,000回計算するのに要した時間を、加算を1flopsとして示した図である。図4において、sqrt(x)は平方根,1/sqrt(x)は逆数平方根をC++言語の標準ライブラリ関数によって演算した場合の計算量を示している。また、fastInvSqrt(x)はニュートン法1回のFast Inverse Square Root法による計算量を示している。fastInvSqrt2(x)はニュートン法2回のFast Inverse Square Root法による計算量を示している。fastInvSqrt3(x)はニュートン法3回のFast Inverse Square Root法による計算量を示している。mFastInvSqrt(x)は修正Fast Inverse Square Root法による計算量を示している。
 また、H.A. Thant, et. al., “Mobile Agents Based Load Balancing Method for Parallel Applications,” 6th Asia-Pacific Symposium on Information and Telecommunication Technologies, Yangon, 2005.によれば、浮動小数点型の加算1回の計算量を同じく1flopとすると、減算、乗算は1flop、除算および平方根計算は4flopsの計算量である。またFastInvSqrt法による逆数平方根演算は、図4に示すとおり、1flopsと見積もることができる。
 a、b、cを従来の計算装置20で計算する場合、計算量の総計は17flopsである。一方、本実施形態に係る計算装置10Aで計算する場合、計算量の総計は10flopsである。よって、従来の計算装置20と比較して、本実施形態に係る計算装置10Aでは、7flops、およそ41.2%の浮動小数点演算を削減すると見積もることができる。
 [第2実施形態]
 以下、第2実施形態における計算装置10である計算装置10Bについて説明する。図5は、実施形態における、和と差の逆数平方根、および和と差の平方根差を計算する計算装置10Bの構成を示すブロック図である。計算装置10Bは、入力装置3から正の実数X、e(X>e>0)が入力され、a=1/√(X+e)、b=1/√(X-e)、c=√(X+e)-√(X-e)を出力装置5に出力する計算装置である。計算装置10Bは、高速逆数平方根計算部31、41、減算部32、乗算部27、37、47、加算部23、33、および除算部44を備える。
 計算装置10Bにおいて、正の実数Xが入力されると、Xは、加算部23および減算部32に入力される。正の実数eが入力されると、eは、加算部23、減算部32、および乗算部27に入力される。加算部23は、X+eを高速逆数平方根計算部31に出力する。減算部32は、X-eを高速逆数平方根計算部41に出力する。乗算部27は、2eを乗算部47に出力する。
 高速逆数平方根計算部31は、乗算部37、加算部33、および出力装置5にa=1/√(X+e)を出力する。高速逆数平方根計算部41は、乗算部37、加算部33、および出力装置5にb=1/√(X-e)を出力する。乗算部37は、ab=(1/√(X+e))(1/√(X-e))を乗算部47に出力する。乗算部47は、2eabを除算部44に出力する。加算部33は、a+b=1/√(X+e)+1/√(X-e)を除算部44に出力する。除算部44は、cとして2eab/(a+b)を出力装置5に出力する。2eab/(a+b)=√(X+e)-√(X-e)であることは、第1実施形態におけるAをX+eとし、BをX-eとすることで示される。
 以上より、計算装置10Bは、a=1/√(X+e)、b=1/√(X-e)、およびc=√(X+e)-√(X-e)を出力装置5に出力する。計算装置10Bは、逆数平方根を計算するが、√(X+e)、√(X-e)は計算していない。よって、cを√(X+e)、√(X-e)を計算することなく出力する。
 図6は、参考例として、従来通り、√(X+e)、√(X-e)を計算してa,b,cを出力する計算装置40を示すブロック図である。計算装置40は、平方根計算部35、45、減算部42、52、加算部43、および除算部54、64を備える。
 計算装置40において、正の実数X、eが入力されると、加算部43は、平方根計算部35に、X+eを出力する。減算部42は、平方根計算部45に、X-eを出力する。平方根計算部35は、除算部54および減算部52に√(X+e)を出力する。平方根計算部45は、除算部64および減算部52に√(X-e)を出力する。
 除算部54は、a=1/√(X+e)を出力する。除算部64は、b=1/√(X-e)を出力する。減算部52は、c=√(X+e)-√(X-e)を出力する。
 a、b、cを従来の計算装置40で計算する場合、計算量の総計は19flopsである。一方、本実施形態に係る計算装置10Bで計算する場合、計算量の総計は12flopsである。よって、従来の計算装置40と比較して、本実施形態に係る計算装置10Bでは、7flops、およそ36.8%の浮動小数点演算を削減すると見積もることができる。
 [第3実施形態]
 以下、第3実施形態における計算装置10である計算装置10Cについて説明する。図7は、実施形態における、M×2行列の第2特異値、特異値の和の逆数および特異値の差の逆数を計算する計算装置10Cの構成を示すブロック図である。本実施形態において、M≧3であり、M=2については第4実施形態で説明する。計算装置10Cは、入力装置3からM×2行列Y=[y,y]が入力され、第2特異値σ、特異値の和の逆数1/(σ+σ)および特異値の差の逆数1/(σ-σ)を出力装置5に出力する計算装置である。
 計算装置10Cは、第2実施形態で説明した計算装置10Bを備える。また、計算装置10Cは、分解部110、内積部16、26、36、減算部62、乗算部57、67、77、87、加算部53、および平方根計算部55を備える。
 計算装置10Cにおいて、M×2行列Y=[y,y]が入力されると、Yは、分解部110に入力される。分解部110は、M×2行列Yを列ベクトルy,yに分解する。分解部110は、yを内積部16、36に出力する。分解部110は、yを内積部26、36に出力する。
 内積部16は、yとyとの内積aを計算し、aを加算部53、乗算部67、および出力装置5に出力する。内積部26は、yとyとの内積cを計算し、cを加算部53、乗算部67、および出力装置5に出力する。内積部36は、yとyとの内積bを計算し、bを乗算部77、および出力装置5に出力する。
 加算部53は、f=a+cを計算装置10Bおよび出力装置5に出力する。乗算部67は、acを減算部62に出力する。乗算部77は、bを減算部62に出力する。減算部62は、d=ac-bを平方根計算部55および出力装置5に出力する。平方根計算部55は、e=√dを乗算部57および出力装置5に出力する。乗算部57は、2eを計算装置10Bに出力する。
 計算装置10Bは、1/√(f+2e)を1/(σ+σ)として出力装置5に出力する。計算装置10Bは、1/√(f-2e)を1/(σ-σ)として出力装置5に出力する。計算装置10Bは、√(f+2e)-√(f-2e)を2σとして乗算部87に出力する。乗算部87は、σを出力装置5に出力する。
 ここで1/√(f±2e)が1/(σ±σ)であることは、以下の通りである。非特許文献1によれば、M×2行列の特異値σ、σの和および差は、σ±σ=√(tr(YY)±2√(det(YY)))である。
 tr(YY)=a(=yとyとの内積)+c(=yとyとの内積)=fである。det(YY)=a(=yとyとの内積)×c(=yとyとの内積)-b(=(yとyとの内積b))=ac-b=dある。よって、√(det(YY))=√d=eである。
 よって、√(tr(YY)±2√(det(YY)))=√(f±2e)である。したがって、1/√(f±2e)=1/(σ±σ)である。
 なお、図7に示されるように、計算装置10Cは、第2特異値σ、特異値の和の逆数1/(σ+σ)および特異値の差の逆数1/(σ-σ)だけではなく、a、b、c、d、e、fも出力装置5に出力する。
 計算装置10Cの計算量について説明する。分解部110の処理は、行列を列ベクトルに分解するだけなので浮動小数点演算は行われない(0flop)。また内積部16、26、36による計算量は、M次列ベクトルの入力に対し、2M-1flopsの計算量である。またsqrtによる平方根計算は前述の通り4flopsである。よって、計算装置10Cによる計算量の総計は6M+19flopsである。上述したように、計算装置10Bにおいて、7flopsの削減効果があるため、計算装置10Cも計算装置10Bによって7flopsの削減効果がある。
 [第4実施形態]
 以下、第4実施形態における計算装置10である計算装置10Dについて説明する。図8は、実施形態における、2×2行列の第2特異値、特異値の和の逆数および特異値の差の逆数を計算する計算装置10Dの構成を示すブロック図である。計算装置10Dは、入力装置3から2×2行列Y=[y,y]が入力され、第2特異値σ、特異値の和の逆数1/(σ+σ)および特異値の差の逆数1/(σ-σ)を出力装置5に出力する計算装置である。
 計算装置10Dは、第2実施形態で説明した計算装置10Bを備える。また、計算装置10Dは、分解部120、内積部46、56、行列式計算部19、絶対値計算部18、乗算部97、107、および加算部63を備える。
 計算装置10Dにおいて、2×2行列Y=[y,y]が入力されると、Yは、分解部120に入力される。分解部120は、2×2行列Yを列ベクトルy,yに分解する。分解部120は、yを内積部46に出力する。分解部120は、yを内積部56に出力する。
 内積部46は、yとyとの内積aを計算し、aを加算部63に出力する。内積部56は、yとyとの内積cを計算し、cを加算部63に出力する。行列式計算部19は、Yの行列式dを計算し、dを絶対値計算部18、および出力装置5に出力する。
 加算部63は、f=a+cを計算装置10Bおよび出力装置5に出力する。絶対値計算部18は、dの絶対値eを乗算部97に出力する。乗算部97は、2eを計算装置10Bに出力する。
 計算装置10Bは、1/√(f+2e)を1/(σ+σ)として出力装置5に出力する。計算装置10Bは、1/√(f-2e)を1/(σ-σ)として出力装置5に出力する。計算装置10Bは、√(f+2e)-√(f-2e)を2σとして乗算部107に出力する。乗算部107は、σを出力装置5に出力する。
 ここで1/√(f±2e)が1/(σ±σ)であることは、第3実施形態で示した通りである。なお、図8に示されるように、計算装置10Dは、第2特異値σ、特異値の和の逆数1/(σ+σ)および特異値の差の逆数1/(σ-σ)だけではなく、d、fも出力装置5に出力する。
 計算装置10Dの計算量について説明する。行列式計算部19の計算量は3flopsである。絶対値計算部18は符号を評価し負の場合に反転させるだけなので浮動小数点演算は行われない(0flop)。計算装置10Dによる計算量の総計は24flopsである。上述したように、計算装置10Bにおいて、7flopsの削減効果があるため、計算装置10Dも計算装置10Bによって7flopsの削減効果がある。
 [第5実施形態]
 以下、第5実施形態における計算装置10である計算装置10Eについて説明する。図9は、実施形態における、M×2行列のSVTを計算する計算装置10Eの構成を示すブロック図である。本実施形態において、M≧3であり、M=2については第6実施形態で説明する。計算装置10Eは、入力装置3からM×2行列Y=[y,y]と正の実数μが入力され、特異値閾値処理(Singular Value Thresholding;SVT)を行い、その特異値閾値処理結果として、下記(1)に示される行列Zを計算する。
Figure JPOXMLDOC01-appb-M000001
 (1)におけるIは2×2の単位行列である。また、γ、δ、e、a、b、c、dについては後述する。
 計算装置10Eは、第3実施形態で説明した計算装置10Cと、係数算出装置200と、M×2行列変換装置300とを備える。
 計算装置10Eにおいて、M×2行列Y=[y,y]が入力されると、Yは、計算装置10C、係数算出装置200、およびM×2行列変換装置300に入力される。また、計算装置10Eにおいて、実数μが入力されると、μは、係数算出装置200に入力される。計算装置10Cは、第3実施形態で説明したように、第2特異値σ、特異値の和の逆数1/(σ+σ)、特異値の差の逆数1/(σ-σ)、a、b、c、d、e、fを出力する。ここで、図9に示されるように、g=1/(σ+σ)、h=1/(σ-σ)とする。
 計算装置10Cの出力のうち、d、f、σ、g、hは、係数算出装置200に出力される。a、b、c、eは、M×2行列変換装置300に出力される。
 係数算出装置200は、d、f、σ、g、h、μが入力される。係数算出装置200は、γ、δを出力する。このγ、δの算出方法について説明する。まず、十分に大きい実数をRとおく。具体的にRの大きさとして、単精度浮動小数点型の最大の数の10分の1程度の大きさが挙げられる。その上で、係数算出装置200は、下記(a)から(d)の4つの場合分けを行うことでγ、δを算出し、それらを出力する。
 (a) Yが零行列の場合
 (b) (a)に該当せず、d=0の場合
 (c) (b)に該当せず、h>Rの場合
 (d) (c)に該当しない場合
 以下、各場合ごとに出力されるγ、δを示す。なお、高速逆数平方根法による逆数平方根をisqrt(・)と表現することがある。例えば、高速逆数平方根法により計算された正の実数Aの逆数平方根を、isqrt(A)と表現することがある。また、(・)における右下添字の+は、ランプ関数を示す。ランプ関数は、入力が負の実数なら0を出力し、入力が非負の実数なら入力された実数をそのまま出力する。関数minは、入力された数値のうちで最も小さい値を出力する。
 (a):γ=0、δ=0
 (b):γ=(1-μ×isqrt(f))、δ=0
 (c):γ=(1-(√2)×μ×isqrt(f))、δ=0
 (d):γ=(1-(μ-σ×h)、δ=min(μ、σ)×g
 係数算出装置200は、上記場合分けにより、γ、δをM×2行列変換装置300に出力する。これにより、M×2行列変換装置300には、上記(1)に含まれるパラメータが全て揃うため、それらを用いて(1)を出力装置5に出力する。なお、e=0の場合には、(1)におけるγδ/eを0とする。
 計算装置10Eの計算量について説明する。ランプ関数および関数minは実数値の比較が主な計算処理であり、浮動小数点演算は行われない(0flop)。図16Bに示されるアルゴリズムに従って計算する場合、計算量は12M+38flopsである。一方、計算装置10Eによる計算量の総計は12M+29flopsである。以上より9flopsの浮動小数点演算を削減すると見積もることができる。
 [第6実施形態]
 以下、第6実施形態における計算装置10である計算装置10Fについて説明する。図10は、実施形態における、2×2行列のSVTを計算する計算装置10Fの構成を示すブロック図である。計算装置10Fは、入力装置3から2×2行列Y=[y,y]と正の実数μが入力され、上述したSVTを行い、その特異値閾値処理結果として、下記(1)に示される行列Zを計算する。なお、yijは、Yのi行j列成分である。
Figure JPOXMLDOC01-appb-M000002
 (2)における関数sign(・)は符号関数であり、入力が負の場合には-1を出力し、入力が0の場合には0を出力し、入力が正の場合に+1を出力する。また、γ、δ、dについては後述する。
 計算装置10Fは、第4実施形態で説明した計算装置10Dと、係数算出装置201と、2×2行列変換装置301とを備える。
 計算装置10Fにおいて、2×2行列Y=[y,y]が入力されると、Yは、計算装置10D、係数算出装置201、およびM×2行列変換装置301に入力される。また、計算装置10Fにおいて、実数μが入力されると、μは、係数算出装置201に入力される。計算装置10Dは、第3実施形態で説明したように、第2特異値σ、特異値の和の逆数1/(σ+σ)、特異値の差の逆数1/(σ-σ)、d、fを出力する。ここで、図10に示されるように、g=1/(σ+σ)、h=1/(σ-σ)とする。
 計算装置10Dの出力のうち、d、f、σ、g、hは、係数算出装置201に出力される。dは、2×2行列変換装置300に出力される。
 係数算出装置201は、d、f、σ、g、h、μが入力される。係数算出装置201は、γ、δを出力する。このγ、δの算出方法について説明する。まず、十分に大きい実数をRとおく。具体的にRの大きさとして、単精度浮動小数点型の最大の数の10分の1程度の大きさが挙げられる。その上で、係数算出装置200は、下記(a)から(d)の4つの場合分けを行うことでγ、δを算出し、それらを出力する。
 (a) Yが零行列の場合
 (b) (a)に該当せず、d=0の場合
 (c) (b)に該当せず、h>Rの場合
 (d) (c)に該当しない場合
 以下、各場合ごとに出力されるγ、δを示す。
 (a):γ=0、δ=0
 (b):γ=(1-μ×isqrt(f))、δ=0
 (c):γ=(1-(√2)×μ×isqrt(f))、δ=0
 (d):γ=(1-(μ-σ×h)、δ=min(μ、σ)×g
 係数算出装置200は、上記場合分けにより、γ、δを2×2行列変換装置301に出力する。これにより、2×2行列変換装置301には、上記(2)に含まれるパラメータが全て揃うため、それらを用いて(2)を出力装置5に出力する。
 計算装置10Fの計算量について説明する。図16Bに示されるアルゴリズムに従って計算する場合、計算量は41flopsである。一方、計算装置10Fによる計算量の総計は32flopsである。以上より9flopsの浮動小数点演算を削減すると見積もることができる。
 次に、グラフの単純化処理について説明する。図11は、グラフの単純化処理の処理概要を示す図である。図11におけるG=(V,E,P)は、M次元空間に埋め込まれたグラフを示す。ここでV,Eはそれぞれ頂点、辺の集合とし、P∈R|V|×Mを頂点座標とする(ここでのRは実数全体の集合)。Pの各行ベクトル(p、(p、…、(p|v|はグラフの各頂点の座標を表す。また、次数が2の頂点の集合をV~={v∈V|degv=2}とする。
 グラフ単純化処理では図11に示される通り、多くの頂点を持つ歪な形状の埋め込みグラフを入力とし、局所的に線形に整列した埋め込みグラフを中間生成し(S101)、最後に不要点を除去して(S102)、形状単純された所望のグラフを得る。
 図12は、グラフ単純化装置500の構成例を示す図である。グラフ単純化装置500は、閾値λとグラフGが入力され、単純化したグラフG’’を出力する。グラフ単純化装置500は、局所線形整列部510と、不要頂点除去部520とを備える。局所線形整列部510は、グラフGを閾値λを用いて局所線形整列させたG’を不要頂点除去部520に出力する。不要頂点除去部520は、入力したグラフG’の頂点のうち、不要な頂点を除去したG’’を出力する。
 図13は、図12における局所線形整列部510の構成例を示す図である。局所線形整列部510は、上述したように、グラフGを閾値λを用いて局所線形整列させたG’を不要頂点除去部520に出力する。局所線形整列部510は、凸最適化問題立式部511と、凸最適化問題求解部512と、座標情報置換部513とを備える。
 最初に凸最適化問題立式部511について説明する。図14は、隣接辺のベクトルを並べる操作を説明するための図である。図14に示される行列Lは、vから隣接する頂点へのベクトル並べるための行列である。図14に示される行列Lの要素について、「…」の箇所は全て0である。また、vの座標がPのk行目としたとき、図14に示されるLにおける「-1」がk列目なっている。LPにより、vを始点としてk-1行目の座標を終点とするベクトルと、vを始点としてk+1行目の座標を終点とするベクトルが得られる。このLを作用させる操作は線形写像である。
 元の入力グラフの形を忠実に再現しながらも曲折回数が少なれば、グラフの局所線形整列化に成功したと言える。ここでは忠実再現の尺度をL1ノルムとし、辺の曲折回数の正則化を核型ノルム関数として、凸最適化問題立式部511は、下記(3)のとおり最適化問題を立式する。ここで核型ノルム関数とは入力行列の特異値の和を計算する関数である。
Figure JPOXMLDOC01-appb-M000003
 上記(3)の最適化問題を解くために、Primal-Dual Splitting(L. Condat, “A primal-dual splitting method for convex optimization involving lipschitzian, proximable and linear composite terms,” Journal of Optimization Theory and Applications, 2013.参照)を用いる。凸最適化問題求解部512が実行する具体的な手順を図15に示す。図15は、局所線形整列問題の解法を示すアルゴリズムを示す図である。なお、図15の4行目のproxτfにおけるfは、下記(4)である。
Figure JPOXMLDOC01-appb-M000004
 図15に示されるアルゴリズムでは、3行目から11行目までのn=1からrまでのr回のループの中に、上述したV~に属する元の全てに対して、6行目から10行目までループが行われることが示されている。
 したがって、6行目から10行目までは、r×(V~に属する元の総数)回実行される。その6行目から10行目のうちの8行目で計算される下記(5)は、核型ノルムの近接写像である。
Figure JPOXMLDOC01-appb-M000005
 上記(5)は、下記(6)の行列の閾値λによるSVTである。
Figure JPOXMLDOC01-appb-M000006
 この8行目の処理は、図15に示されるアルゴリズムの計算時間のうちの約53.8%を占める。そこで、上述した計算装置10E、10Fを用いて上記SVTを計算することで、従来と比較して、グラフ単純化処理を高速に実行することができる。
 座標情報置換部513は、凸最適化問題求解部512により局所線形整列された頂点座標に座標情報を置換して、局所線形整列させたG’を出力する。
 本実施形態は、グラフの単純化だけではなく、SVTを行う全ての処理に適用可能である。例えば、画像偽色除去は、グラフの単純化と同様に、多数の小型行列を正則化する問題に分類される。
 以上説明したように、本実施形態によれば、核型ノルムを特異値を用いずに混合ノルムで表現することで、SVDが不要なSVT計算を実現することで、計算量を削減可能となる。さらに、アルゴリズムを容易にデータ並列化でき、多数の行列について同時処理が可能である。アルゴリズムをデータ並列化できれば、パソコンに搭載されるCPUの多くが採用しているSingle Instrucion Multiple Data(SIMD)等のデータ並列アーキテクチャを用いる実装で高速化できる。
 以上説明した計算装置10A~10Fにおいて、加算部や内積部などの各種計算を行う構成が計算結果を一時的に記憶するメモリなどの記憶装置を設け、この記憶装置に計算結果を一時的に記憶してもよい。
 以上、この発明の実施形態について図面を参照して詳述してきたが、具体的な構成はこの実施形態に限られるものではなく、この発明の要旨を逸脱しない範囲の設計等も含まれる。
 本発明は、平方根の計算を行う計算装置に適用可能である。
10、10A、10B、10C、10D、10E、10F、20、40…計算装置、11、21、31、41…高速逆数平方根計算部、12、22、32、42、52…減算部、13、23、33、43、53、63…加算部、14、24、34、44、54、64…除算部、15、25、35、45、55…平方根計算部、16、26、36、46、56…内積部、17、27、37、47、57、67、77、87、97、107、117…乗算部、18…絶対値計算部、19…行列式計算部、110、120…分解部、200、201…係数算出装置、300、301…行列変換装置、500…グラフ単純化装置、510…局所線形整列部、511…凸最適化問題立式部、512…凸最適化問題求解部、513…座標情報置換部、520…不要頂点除去部

Claims (7)

  1.  コンピュータが、正の実数Aの逆数平方根を高速逆数平方根法により計算する第1逆数平方根ステップと、
     コンピュータが、正の実数Bの逆数平方根を高速逆数平方根法により計算する第2逆数平方根ステップと、
     コンピュータが、AからBを減算する減算ステップと、
     コンピュータが、前記第1逆数平方根ステップでの計算結果と、前記第2逆数平方根ステップでの計算結果とを乗算する第1乗算ステップと、
     コンピュータが、前記第1乗算ステップでの計算結果と、前記減算ステップでの計算結果とを乗算する第2乗算ステップと、
     コンピュータが、前記第1逆数平方根ステップでの計算結果と、前記第2逆数平方根ステップでの計算結果とを加算する加算ステップと、
     コンピュータが、前記第2乗算ステップでの計算結果を、前記加算ステップでの計算結果で除算する除算ステップと、
     を備えた計算方法。
  2.  コンピュータが、正の実数Xと当該実数Xよりも小さい正の実数eとを加算する第1加算ステップと、
     コンピュータが、Xからeを減算する減算ステップと、
     コンピュータが、eを2倍する第1乗算ステップと、
     コンピュータが、前記第1加算ステップでの計算結果の逆数平方根を高速逆数平方根法により計算する第1逆数平方根ステップと、
     コンピュータが、前記減算ステップでの計算結果の逆数平方根を高速逆数平方根法により計算する第2逆数平方根ステップと、
     コンピュータが、前記第1逆数平方根ステップでの計算結果と、前記第2逆数平方根ステップでの計算結果とを加算する第2加算ステップと、
     コンピュータが、前記第1逆数平方根ステップでの計算結果と、前記第2逆数平方根ステップでの計算結果とを乗算する第2乗算ステップと、
     コンピュータが、前記第2乗算ステップでの計算結果と、前記減算ステップでの計算結果とを乗算する第3乗算ステップと、
     コンピュータが、前記第3乗算ステップでの計算結果を、前記第2加算ステップでの計算結果で除算する除算ステップと、
     を備えた計算方法。
  3.  コンピュータが、M×2行列を構成する第1列ベクトルと第2列ベクトルのうち、第1列ベクトル同士の内積を計算する第1内積ステップと、
     コンピュータが、前記第2列ベクトル同士の内積を計算する第2内積ステップと、
     コンピュータが、第1列ベクトルと第2列ベクトルとの内積を計算する第3内積ステップと、
     コンピュータが、前記第1内積ステップでの計算結果と、前記第2内積ステップでの計算結果とを加算する加算ステップと、
     コンピュータが、前記第1内積ステップでの計算結果と、前記第2内積ステップでの計算結果とを乗算する第1乗算ステップと、
     コンピュータが、前記第3内積ステップで計算された内積同士を乗算する第2乗算ステップと、
     コンピュータが、前記第1乗算ステップでの計算結果から前記第2乗算ステップでの計算結果を減算する減算ステップと、
     コンピュータが、前記減算ステップでの計算結果の平方根を計算する平方根計算ステップと、
     コンピュータが、前記平方根計算ステップでの計算結果を2倍する第3乗算ステップと、
     を備え、
     前記加算ステップでの計算結果をXとし、第3乗算ステップでの計算結果をeとして、請求項2に記載の計算方法により、M×2行列の2つの特異値の和の逆数、2つの特異値の差の逆数、および1つの特異値を計算する計算方法。
  4.  請求項3の計算方法より計算された2つの特異値の和の逆数、および2つの特異値の差の逆数を用いて、コンピュータが、M×2行列の特異値閾値処理を計算するステップを備えた計算方法。
  5.  コンピュータが、M次元空間のグラフの頂点を局所的に線形に整列する整列ステップと、
     コンピュータが、前記整列ステップにおいて線形に整列されたグラフの頂点のうち、頂点同士を結ぶ辺のなす角度に基づき、不要な辺と頂点とを除去する除去ステップと、
     を備え、
     前記整列ステップにおいて、特異値の和を計算する核型ノルム関数による特異値閾値処理を前記請求項4に記載の計算方法で行う計算方法。
  6.  請求項1から請求項5のいずれか1項に記載の計算方法をコンピュータに実行させるためのプログラム。
  7.  正の実数Aの逆数平方根を高速逆数平方根法により計算する第1逆数平方根部と、
     正の実数Bの逆数平方根を高速逆数平方根法により計算する第2逆数平方根部と、
     AからBを減算する減算部と、
     前記第1逆数平方根部での計算結果と、前記第2逆数平方根部での計算結果とを乗算する第1乗算部と、
     前記第1乗算部での計算結果と、前記減算部での計算結果とを乗算する第2乗算部と、
     前記第1逆数平方根部での計算結果と、前記第2逆数平方根部での計算結果とを加算する加算部と、
     前記第2乗算部での計算結果を、前記加算部での計算結果で除算する除算部と、
     を備えた計算装置。
PCT/JP2020/043722 2020-11-25 2020-11-25 計算方法、計算装置およびプログラム WO2022113180A1 (ja)

Priority Applications (2)

Application Number Priority Date Filing Date Title
PCT/JP2020/043722 WO2022113180A1 (ja) 2020-11-25 2020-11-25 計算方法、計算装置およびプログラム
JP2022564862A JP7560765B2 (ja) 2020-11-25 2020-11-25 計算方法、計算装置およびプログラム

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
PCT/JP2020/043722 WO2022113180A1 (ja) 2020-11-25 2020-11-25 計算方法、計算装置およびプログラム

Publications (1)

Publication Number Publication Date
WO2022113180A1 true WO2022113180A1 (ja) 2022-06-02

Family

ID=81754217

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2020/043722 WO2022113180A1 (ja) 2020-11-25 2020-11-25 計算方法、計算装置およびプログラム

Country Status (2)

Country Link
JP (1) JP7560765B2 (ja)
WO (1) WO2022113180A1 (ja)

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20020039454A1 (en) * 2000-06-15 2002-04-04 Lifef/X Networks, Inc. And Auckland Uniservices Li Basis functions of three-dimensional models for compression, transformation and streaming
JP2019046196A (ja) * 2017-09-01 2019-03-22 日本電信電話株式会社 行列単純化装置、プログラム、および行列単純化方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20020039454A1 (en) * 2000-06-15 2002-04-04 Lifef/X Networks, Inc. And Auckland Uniservices Li Basis functions of three-dimensional models for compression, transformation and streaming
JP2019046196A (ja) * 2017-09-01 2019-03-22 日本電信電話株式会社 行列単純化装置、プログラム、および行列単純化方法

Also Published As

Publication number Publication date
JPWO2022113180A1 (ja) 2022-06-02
JP7560765B2 (ja) 2024-10-03

Similar Documents

Publication Publication Date Title
US11093669B2 (en) Method and system for quantum computing
Bauer et al. Clear and compress: Computing persistent homology in chunks
Özbay et al. Poisson CNN: Convolutional neural networks for the solution of the Poisson equation on a Cartesian mesh
WO2019168084A1 (ja) 推論装置、畳み込み演算実行方法及びプログラム
WO2015125714A1 (en) Method for solving convex quadratic program for convex set
EP3507747A1 (en) Exact quantum circuits and circuit syntheses for qudit and multiple qubit circuits
US20190347072A1 (en) Block floating point computations using shared exponents
JP2015197702A (ja) 情報処理装置、情報処理方法
WO2020221583A1 (en) System and method for molecular design on a quantum computer
Diviánszky et al. Qutrit witness from the Grothendieck constant of order four
CN112634149A (zh) 一种基于图卷积网络的点云去噪方法
Koehl et al. Statistical physics approach to the optimal transport problem
Peng et al. Mbfquant: a multiplier-bitwidth-fixed, mixed-precision quantization method for mobile cnn-based applications
WO2020223850A1 (en) System and method for quantum circuit simulation
JP6810003B2 (ja) 行列単純化装置、プログラム、および行列単純化方法
Bizzarri et al. Symmetries of discrete curves and point clouds via trigonometric interpolation
Zhou et al. Quantum multidimensional color images similarity comparison
WO2022113180A1 (ja) 計算方法、計算装置およびプログラム
US11631002B2 (en) Information processing device and information processing method
JP7091930B2 (ja) テンソルデータ計算装置、テンソルデータ計算方法及びプログラム
Doukhnitch et al. Hardware-oriented algorithm for quaternion-valued matrix decomposition
WO2024027933A1 (en) Quantum computing method for solving combinatorial optimization problems
CN113924610B (zh) 秘密共轭梯度法计算系统及方法、秘密计算装置、共轭梯度法计算装置及方法、以及记录介质
Ramadas et al. Testing nelder-mead based repulsion algorithms for multiple roots of nonlinear systems via a two-level factorial design of experiments
WO2020040007A1 (ja) 学習装置、学習方法及び学習プログラム

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

Country of ref document: EP

Kind code of ref document: A1

ENP Entry into the national phase

Ref document number: 2022564862

Country of ref document: JP

Kind code of ref document: A

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 20963445

Country of ref document: EP

Kind code of ref document: A1