It is common to calculate sum. Usually, a sum over $n$ numbers has $n - 1$ roundings.
Use $fl(x)$ as the floating number round $x$ (usually rounded by hardware, the closest representable number, with some even/odd considerations). So $n$ numbers $x_1, x_2, \ldots, x_n$ is usually calculated as
$$S(x_1, x_2, \ldots, x_n) = fl(x_n + fl(x_{n-1} + \cdots ))$$
This calculation is acceptable in most case since floating type, like double, has quite large precision. Things become worse when the number of values are large, and values are in a wild range, esp. when a lot of them could cancel each other, and the result is very small compared to some numbers. In tradition, we could use Kahan summation algorithm to improve the precision.
In this post, I try a simple algorithm which would give the rounding floating point number of the exact sum.
Two number sum could be represented exactly as the floating sum and residual, i.e. $a + b = fl(a + b) + e$, with e also a floating number. In Knuth's TAOCP book, there is a simple algorithm discussion to find the $e$.
Similar to this idea, and also considering floating formats usually have fixed width digits (mantissa part) and exponents, so the sum could be represented by a limit numbers sum. That is, for numbers $x_1, x_2, \ldots, x_n$, we could find up to $L$ numbers $y_j$, such that
$$\sum_i x_i = \sum_{j=1}^{L} y_j$$.
We can assume that a minimal number of $y_i$ have no zeros,and no overlapped precession, i.e. any $i$, $j$, $fl(y_i + y_j) = y_i$, or $y_j$, since we could always rewrite $y_i + y_j$ to $fl(y_i + y_j) + e$. Based on this, we could bound the $L$ as the exponent span divided by the digits size.
To calculate numbers sum, we manage a number series (length up to $L$), starting from empty. We add numbers to this series one by one. When the length exceeds the limit $L$, we merge some overlaps $y_i + y_j$ to $fl(y_i + y_j) + e$, and some $e$ or even $fl(y_i + y_j)$ are zeros and omitted, which makes the length bounded by $L$ again. After all numbers added, we merge all overlaps and the number in the series with largest absolute value is the final result.
Code:
Download all pdf files from Google scholar alert emails
I have gathered two many Google scholar alert emails. When i am with my laptops, i have much more fun things to do than reading papers. But when I am with my e-reader, there are only a small number of older documents. To keep documents in the e-reader updated, i added a cronjob to parse the scholar alert emails in Gmail's inbox, download all pdf files, push them into Google drive, so that next time, i could simply sync these folders again with my e-reader.
The alert email parsing file is a pure Python script, but the shell script which downloads gmails, and syncs files between local and Google drive depending on some personal binaries, which need to be replaced if you want to use, sorry. The csv_sql could be simply replaced by sqlite3, the pjobs could be replaced by parallel from moreutils, and sync between local and Google drive is optional.
The alert email parsing file is a pure Python script, but the shell script which downloads gmails, and syncs files between local and Google drive depending on some personal binaries, which need to be replaced if you want to use, sorry. The csv_sql could be simply replaced by sqlite3, the pjobs could be replaced by parallel from moreutils, and sync between local and Google drive is optional.
A naive secure model of secure sharing
I have used secure sharing to distribute my private encrypted data for many years. And i have a demo project in my github. FYI, this demo project is not the one used by me in these years, since it is just a demo and is not safe.
In a secure sharing, we say every sharing as a part of the original message, and we can denote the $i$-th part by $P_i$.
To protect our data, we have to prove the secure sharing scheme is safe. A naive definition of sharing safe can be:
The above text can transferred to, for any given message distribution, and some parts $P_{i_j}$ of any message extracted from this distribution,
$$P(P_{i_0}|P_{i_1}, P_{i_2}, \ldots, P_{i_{K-1}}) = P(P_{i_0}).$$
A easiest way to achieve this is using a weaker secure shared transformation T, which is not part of message, and based on the transformation, we have, any (K - 1) parts are valid for any message. That is, for any message M, and any (K - 1) parts $P_0, P_1, \ldots, P_{K-1}$, we have a transformation $T$, such that $S_i(M|T) = P_i$, where $S_i$ is a $i$-th sharing part under transformation $T$. If we can get the $T$ from $K - 1$ parts, the secure is still not guaranteed.
But it's easy to prove that a revertible transformation can be secured shared in a sense that:
That is, we lost one random dimension. If this is not acceptable, we can have a chain of transformations, and this chain will converge to the real random secure model.
If you have different views of this secure sharing model, please kindly let me know, so that i'm not in a risk I do not know.
Terms
In a secure sharing, we say every sharing as a part of the original message, and we can denote the $i$-th part by $P_i$.
Probability Model
To protect our data, we have to prove the secure sharing scheme is safe. A naive definition of sharing safe can be:
Given any message distribution, the likelihood of a part is independent of any other (K - 1) parts.
The above text can transferred to, for any given message distribution, and some parts $P_{i_j}$ of any message extracted from this distribution,
$$P(P_{i_0}|P_{i_1}, P_{i_2}, \ldots, P_{i_{K-1}}) = P(P_{i_0}).$$
A easiest way to achieve this is using a weaker secure shared transformation T, which is not part of message, and based on the transformation, we have, any (K - 1) parts are valid for any message. That is, for any message M, and any (K - 1) parts $P_0, P_1, \ldots, P_{K-1}$, we have a transformation $T$, such that $S_i(M|T) = P_i$, where $S_i$ is a $i$-th sharing part under transformation $T$. If we can get the $T$ from $K - 1$ parts, the secure is still not guaranteed.
But it's easy to prove that a revertible transformation can be secured shared in a sense that:
Given any $K - 1$ rows of a revertible transformation matrix, the space for the last row is isomorphic to $F^{K-1}$, where $F$ is the under field space.
That is, we lost one random dimension. If this is not acceptable, we can have a chain of transformations, and this chain will converge to the real random secure model.
Practice
In practice, only one transformation, plus a random accumulated random vector, give quite high entropy of every parts, which is verified by gzip.If you have different views of this secure sharing model, please kindly let me know, so that i'm not in a risk I do not know.
CMake to mimic Google Build language
Google use a description language for build and dependency management. That is very convenient, in order to use some libraries, you just need to include some headers in your C++, and add that libraries path to 'deps' list.
I'm using CMake for personal tiny projects. One CMake feature that annoys me a lot is no target namespaces. We can not define the same target in different directories. So I thought, we can define some Google BUILD like functions for CMake and make every target dependent on the it's directory. In my personal tiny projects, I'm using weizi as root of source directory, and in every sub directory, we can refer targets in current directory directly and targets in other directory with full path.
We can not use '/' in target names directly, one choice is using '.' to join directories.
I'm using CMake for personal tiny projects. One CMake feature that annoys me a lot is no target namespaces. We can not define the same target in different directories. So I thought, we can define some Google BUILD like functions for CMake and make every target dependent on the it's directory. In my personal tiny projects, I'm using weizi as root of source directory, and in every sub directory, we can refer targets in current directory directly and targets in other directory with full path.
CC_BINARY(diff_spell_check
SRCS diff_spell_check.cc
LIBS /options_parser/options_parser)
CC_BINARY(csv2sqlite
SRCS csv2sqlite.cc
LIBS /options_parser/options_parser
/third_party/sqlite)
CC_TEST(csv_split_test
SRCS csv_split_test.cc
LIBS /test/test_lite_main)
We can not use '/' in target names directly, one choice is using '.' to join directories.
Real symmetric matrix exponent
General matrix exponent is not trivial, scipy's implementation(help scipy.linalg.expm) refers a paper A New Scaling and Squaring Algorithm for the Matrix Exponential, and here is a revisited version.
But in most optimization cases, the matrices are real and symmetric, even positive definite. We can build a easy to understand method based on SVD. Let us assume that
$$A = USV'$$
is the SVD decomposition of a symmetric matrix A, we have
$$A = USV^T, A^T = VSU^T = A$$
$$A^2 = USV^TVSU^T = US^2U^T$$
$$A^3 = US^2U^TUSV^T = US^3V^T$$
$$A^{2k} = US^{2k}U^T$$
$$A^{2k+1} = US^{2k+1}V^T$$
This suggests that
$$e^A = \sum_k \dfrac{A^k}{k!} = \sum_k U\dfrac{S^{2k}}{(2k)!}U^T + \sum_k U\dfrac{S^{2k+1}}{(2k+1)!}V^T$$
Which leads to
$$e^A = \dfrac{U(e^S + e^{-S})U^T + U(e^S - e^{-S})V^T}{2}$$
If A is a positive definite matrix, then $U=V$, so the above formula can be simplified as $e^A=Ue^SU^T$.
But in most optimization cases, the matrices are real and symmetric, even positive definite. We can build a easy to understand method based on SVD. Let us assume that
$$A = USV'$$
is the SVD decomposition of a symmetric matrix A, we have
$$A = USV^T, A^T = VSU^T = A$$
$$A^2 = USV^TVSU^T = US^2U^T$$
$$A^3 = US^2U^TUSV^T = US^3V^T$$
$$A^{2k} = US^{2k}U^T$$
$$A^{2k+1} = US^{2k+1}V^T$$
This suggests that
$$e^A = \sum_k \dfrac{A^k}{k!} = \sum_k U\dfrac{S^{2k}}{(2k)!}U^T + \sum_k U\dfrac{S^{2k+1}}{(2k+1)!}V^T$$
Which leads to
$$e^A = \dfrac{U(e^S + e^{-S})U^T + U(e^S - e^{-S})V^T}{2}$$
If A is a positive definite matrix, then $U=V$, so the above formula can be simplified as $e^A=Ue^SU^T$.
Safer/slower tape backup
Last week, I read an nice article relative to a video talking about Google backup: How Google Backs Up the Internet.
The article tells us that one way to back up data is using tapes. For every 4 tapes, a fifth tape, the xor of the 4 tapes, is used to recover from tapes corruption.
This strategy can recovers from any one tape corruption of 5 tapes, and about 20% capacity of tape is wasted. If we assumes simply that every tape have the same corrupt probability p in a fixed period(The time to next check), then if we have 100 tapes, the probability without data lost($P_1$) is $(4p(1-p)^4 + (1-p)^5)^{20}$.
Here We can assume that every tape is self checked, so we can check tapes independently.
The mainly problem is that, if there are 2 tapes corruption in one group of 5 tapes, we'll lost data. So we still have a big chance to lost some data.
I was thinking, there should be a another better strategy. We just need a encoding way, for a given number of tapes(n), use some extra number of tapes(e), so that we can recover the data at a limit number(k, $k \le e$) of tapes lost. The efficiency is $\frac{n}{n + e}$, and the probability without data lost($P_2$) is $\sum_{i=0}^{k}{n + e \choose i}(1-p)^{n + e - i}p^i$.
Or, let $X = 1$ if the tape is corrupted, otherwise $X=0$, and $S = \frac{\sum X}{n + e}$, The above can be rewrite as $S \le \frac{k}{n + e}$. Since $E\; S = E X = p$, $\mathrm{var}\;S = \frac{p(1-p)}{n + e}$, and $S \to p$. If we can keep $\frac{k}{n + e}$, and assume $\frac{k}{n + e} \gt p$, then the larger $n$, $e$, the safer we are. And for large $n, e$, we can assume $S$ suits a normal distrubution, and
\begin{align*}
P\left(S \le \frac{k}{n + e}\right) & = \Phi\left(\frac{\frac{k}{n + e} - \mu}{\sigma}\right) \\
& = \Phi\left(\frac{\sqrt{n + e}(\frac{k}{n + e} - p)}{\sqrt{p (1 - p)}}\right) \\
& \to 1
\end{align*}
Usually, $p$ is very small. For example, if $p=10^{-5}$, then for 100 tapes, 80 data tapes and 20 redundant tapes, we'll survival in $P_1 = 0.9998$. And if we use $n=80$, $k=e=20$, we'll survival in $P_2 = 1.0 - 2.04\times 10^{-84}$.
Keep the same tapes configuration, for different $p$, we get this numbers:
The strategy for this safer way does exist, for example, polynomial code, almost the reverse of secret sharing. We build the virtual whole message from first $n$ data parts, and split the whole message to $n + e$ parts, the first $n$ parts should be identity to original $n$ data parts. The problem is that, in order to recover a tape, we need at least $n$ tapes.
I once heard that, disk vendors have very good redundant encoding algorithms, what are they? May this problem have been resolved very well in that domain, we just need to transfer the solutions.
The article tells us that one way to back up data is using tapes. For every 4 tapes, a fifth tape, the xor of the 4 tapes, is used to recover from tapes corruption.
This strategy can recovers from any one tape corruption of 5 tapes, and about 20% capacity of tape is wasted. If we assumes simply that every tape have the same corrupt probability p in a fixed period(The time to next check), then if we have 100 tapes, the probability without data lost($P_1$) is $(4p(1-p)^4 + (1-p)^5)^{20}$.
Here We can assume that every tape is self checked, so we can check tapes independently.
The mainly problem is that, if there are 2 tapes corruption in one group of 5 tapes, we'll lost data. So we still have a big chance to lost some data.
I was thinking, there should be a another better strategy. We just need a encoding way, for a given number of tapes(n), use some extra number of tapes(e), so that we can recover the data at a limit number(k, $k \le e$) of tapes lost. The efficiency is $\frac{n}{n + e}$, and the probability without data lost($P_2$) is $\sum_{i=0}^{k}{n + e \choose i}(1-p)^{n + e - i}p^i$.
Or, let $X = 1$ if the tape is corrupted, otherwise $X=0$, and $S = \frac{\sum X}{n + e}$, The above can be rewrite as $S \le \frac{k}{n + e}$. Since $E\; S = E X = p$, $\mathrm{var}\;S = \frac{p(1-p)}{n + e}$, and $S \to p$. If we can keep $\frac{k}{n + e}$, and assume $\frac{k}{n + e} \gt p$, then the larger $n$, $e$, the safer we are. And for large $n, e$, we can assume $S$ suits a normal distrubution, and
\begin{align*}
P\left(S \le \frac{k}{n + e}\right) & = \Phi\left(\frac{\frac{k}{n + e} - \mu}{\sigma}\right) \\
& = \Phi\left(\frac{\sqrt{n + e}(\frac{k}{n + e} - p)}{\sqrt{p (1 - p)}}\right) \\
& \to 1
\end{align*}
Keep the same tapes configuration, for different $p$, we get this numbers:
| $p$ | $P_1$ | $P_2$ |
|---|---|---|
| 0.2 | 0.0002 | 0.559 |
| 0.1 | 0.042 | 0.9992 |
| 0.01 | 0.81 | $1.0 - 9.58\times 10^{-22}$ |
| 0.001 | 0.98 | $1.0 - 1.9\times 10^{-42}$ |
| 0.0001 | 0.998 | $1.0 - 2.03\times 10^{-63}$ |
| 0.00001 | 0.9998 | $1.0 - 2.04\times 10^{-84}$ |
| 0.000001 | 0.99998 | $1.0 - 2.04\times 10^{-105}$ |
The strategy for this safer way does exist, for example, polynomial code, almost the reverse of secret sharing. We build the virtual whole message from first $n$ data parts, and split the whole message to $n + e$ parts, the first $n$ parts should be identity to original $n$ data parts. The problem is that, in order to recover a tape, we need at least $n$ tapes.
I once heard that, disk vendors have very good redundant encoding algorithms, what are they? May this problem have been resolved very well in that domain, we just need to transfer the solutions.
Optimization Approach for Boosting
In the Adaptive Boosting paper, the authors gave a construction way to get a stronger learner from weak learners.
If we free the construction from AdaBoost, we can consider the AdaBoost as optimization over linear function spaces spanned by weak learners, the construction method is like a coordinate decent solution algorithm.
We try to give some theory limits of the optimal point under some reasonable(?) assumptions.
Let $w_i$, $i=0, 1, \ldots, n$ be weaker learners. They are totally independent. For any sample $x$, let $y(x)$ be the real class label(-1, 1, or real value). The sign of $w_i(x)y(x)$ will indicate the correctness of $w_i$ on $x$. In most case, $w_i(x)$ has this form $w_i(x) = \sum_j^\infty f_{i_j} \times x_{i_j}$, where $x_{i_j}$ are some features of $x$, transferred by an unknown kernel function of $w_i$.
Here, we assume that the value of $w_i(x) y(x)$ follows $N(\mu_i, \sigma_i^2)$. Then $\Phi(\frac{\mu_i}{\sigma_i})$ is the precision of $w_i$. We also assume that all $w_i(x)y(x)$ are independent, then $\sum_i\alpha_i w_i(x) y(x)$ follows $N(\sum_i \alpha_i\mu_i, \sum_i \alpha_i^2\sigma_i^2)$.
The best $\alpha_i$ is $\mathrm{arg max}_{\alpha}\Phi(\frac{\sum_i\alpha_i\mu_i}{\sqrt{\sum_i\alpha_i^2\sigma_i^2}})$, which can be calculated as $\alpha_i = \frac{c \mu_i}{\sigma_i^2}$, where $c$ is a non-zero constant.
If we assume that, all weak learners has the same precision $p$, we can calculate the stronger learner has precision $\Phi(\sqrt{n}\Phi^{-1}(p))$.
I think this approach is easier to understand than original paper and easier for me to implement with optimization framework. We also has a good start point. We can try different lost and link functions for boosting. But this approach may introduced more over-fitting compared to construction method of original Ada Boosting, which it's not analyzed in this article.
Experiment
For easy of experiment, we draw $n$ samples from $U(0, 1)^d$, and let real positive label as $\mathrm{mean}(x) > 0.5$, we choose some weak learners as $x_i > 0.7$, and $\mathrm{mean}(x) > 0.45$. Note that these weak learners are not independent, but we still got the best weights by the calculation. There're no precision differences between lost functions $e^{-w(x)y(x)}$, $\log(1 + e^{-w(x)y(x)})$ and our direct calculation. You can check code in boosting.go, and run it as following:git clone https://github.com/zuoyan/optimization_go cd optimization_go export GOPATH=$(pwd) go run bin/boosting.go --generate_samples=10000 --generate_dimension=10 --generate_output=train.labels go run bin/boosting.go --generate_samples=10000 --generate_dimension=10 --generate_output=test.labels --seed=13 go run bin/boosting.go --train_file=train.labels --test_file=test.labels --calculate_weights --optimize_weights --eval_weakers
And once we got this(this can regenerate, since math.rand in golang is seeded):
weak precisions at train 0.5968 0.6022 0.5954 0.5947 0.5959 0.5995 0.5929 0.5872 0.6011 0.5995 0.7933 weak precisions at test 0.5918 0.5895 0.5869 0.5975 0.5867 0.5948 0.5931 0.597 0.59 0.5922 0.7953 weak precisions at all 0.5943 0.59585 0.59115 0.5961 0.5913 0.59715 0.593 0.5921 0.59555 0.59585 0.7943 calculate weights 0.06936016186937585 0.07355792072038556 0.06828066537841014 0.06774223410516975 0.06866579457265785 0.07145217471411021 0.06636165980817825 0.06202615041687599 0.07269833592484762 0.07145217471411021 0.30840272777587846 weights precision at train 0.834 weights precision at test 0.8357 weights precision at all 0.83485 2014/04/05 08:16:54 optimize ... 2014/04/05 08:16:54 solver[level=100]:iter=0 y=0.6207724781251298 #f=2/3 #g=2/3 result=Forward 2014/04/05 08:16:54 solver[level=100]:iter=1 y=0.5506796109437893 #f=1/4 #g=1/4 result=Forward 2014/04/05 08:16:54 solver[level=100]:iter=2 y=0.533068799117156 #f=2/6 #g=2/6 result=Forward 2014/04/05 08:16:54 solver[level=100]:iter=3 y=0.5319598811608891 #f=2/8 #g=2/8 result=Forward 2014/04/05 08:16:54 solver[level=100]:iter=4 y=0.531943122354469 #f=2/10 #g=2/10 result=Forward 2014/04/05 08:16:54 solver[level=100]:iter=5 y=0.5319426989350947 #f=2/12 #g=2/12 result=Forward 2014/04/05 08:16:54 solver[level=100]:iter=6 y=0.5319426852247808 #f=2/14 #g=2/14 result=Forward 2014/04/05 08:16:54 solver[level=100]:iter=7 y=0.5319426852122358 #f=2/16 #g=2/16 result=Forward 2014/04/05 08:16:54 solver[level=100]:iter=8 y=0.5319426852121972 #f=1/17 #g=1/17 result=Forward 2014/04/05 08:16:54 solver[level=100]:iter=9 y=0.5319426852121972 #f=3/20 #g=3/20 result=BreakKeep 2014/04/05 08:16:54 solver min value 0.5319426852121972 #f=20 #g=20 at [0.3588231913006707, 0.3626950627657966, 0.3320246202662729, 0.3396292934695077, 0.32662572027393794, 0.32597136863617504, 0.32435762645257327, 0.3182808901882726, 0.3276780489197852, 0.34074165915815785, 1.67455936635565] optimize weights [0.07131695537553615 0.07208649895909681 0.06599067619146069 0.06750212292241123 0.06491763208738632 0.0647875781564208 0.06446684309222839 0.06325907743076396 0.06512678488713779 0.067723208225987 0.33282262267157076] weights precision at train 0.834 weights precision at test 0.8357 weights precision at all 0.83485
Subscribe to:
Posts (Atom)