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.

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$.

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:


$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

Mimic the interface of golang -- Part 1

I'm learning golang. As a C/C++ programmer, I like the feature of open method and dynamic interface, which are not supported(or not supported in native) by C/C++. In my experience,  dynamic interface damage performance a lot, so i write a small bench code:

package main
import (
"fmt"
"time"
)
var (
count = 0
)
func benchmarkFunc(N int, f func(int)) float64 {
start := time.Now()
f(N)
return float64(time.Since(start).Nanoseconds()) / float64(N)
}
func benchmarkFuncPrint(name string, N int, f func(int)) {
ns := benchmarkFunc(N, f)
fmt.Printf("%s %f nano-seconds/call\n", name, ns)
fmt.Printf("%s %f #/second\n", name, 1.e9/ns)
}
type Showable interface {
Show()
}
type Show1 struct {
}
func (p *Show1) Show() {
count++
}
type Show2 struct {
}
func (p *Show2) Show() {
count += 2
}
func call(p interface{}) {
p.(Showable).Show()
}
func call_show(s Showable) {
s.Show()
}
func dynamic(p interface{}, N int) {
for i := 0; i < N; i++ {
// call(p)
p.(Showable).Show()
}
}
func static(s Showable, N int) {
for i := 0; i < N; i++ {
// call_show(s)
s.Show()
}
}
func main() {
N := 10000000
benchmarkFuncPrint("dynamic", N,
func(N int) {
d := Show1{}
dynamic(&d, N)
})
benchmarkFuncPrint("static", N,
func(N int) {
d := Show2{}
static(&d, N)
})
}
This gives
dynamic 19.478761 nano-seconds/call
dynamic 51337966.752537 #/second
static 2.507811 nano-seconds/call
static 398754196.190095 #/second
Which is very quick. Only 19 nano seconds, in my last benchmark, only mutex will cost about 10 nano seconds. And the static version is almost the same as virtual function in C++.

I wrote a toy version of open method for C++, using std::unorderd_map, std::map, google::dense_hash_map, as the storage from std::type_info* to function pointer. And this's the performance counters
lambda: 0.795678 nano-seconds/call
lambda: 1.25679e+09 #/second
unordered_map: 30.4666 nano-seconds/call
unordered_map: 3.28228e+07 #/second
dense hash map: 22.4284 nano-seconds/call
dense hash map: 4.45863e+07 #/second
map: 22.5974 nano-seconds/call
map: 4.42528e+07 #/second
std::map is quicker in this small test, but in larger project, we should divide the performance counter by log(N), where N is the number of function pointers.

It seems that golang's runtime is quicker than google::dense_hash_map + std::mutex. So I felt interesting about it's implementation, and checked the file src/pkg/runtime/iface.c. golang uses a static size(1009) hash table, and only lock if finding first time failed. Chain list in every bucket of hash table is a multiple readers and single writer forward list, reader is wait free.

It's a good implementation, but I don't like the fix size to 1009?(consider a larger project, your performance counter should be divided by ceil(N / 1009), where N is the number of interfaces. And, I guess, using a bit modulo with quadratic open addressing may improve the performance.

So I tried, and I think it did(only test in C++).
template
struct ValueFree {
  inline void operator()(const T & v) const {}
};
template <>
struct ValueFree {
  inline void operator()(void (*v)(void)) const {  }
};

template
struct ValueFree {
  inline void operator()(T * v) const {
    delete v;
  }
};
template
struct TableHash : std::hash {};
template
struct TableHash {
  inline size_t operator()(const T*ptr) const {
    return (intptr_t)(ptr) >> 4;
  }
};
template ,
          class DeleteValue=ValueFree>
struct Table {
  struct State {
    State * old;
    size_t size;
    size_t buckets;
    std::pair table[0];
  };
  std::atomic state_;
  std::mutex add_mutex_;
  std::atomic num_find_;
  Table() {
    const size_t buckets = 8;
    auto state = (State*)malloc(sizeof(State) + sizeof(state_.load()->table[0]) * buckets);
    state->old = nullptr;
    state->size = 0;
    state->buckets = buckets;
    state_.store(state);
    num_find_.store(0);
  }
  ~Table() {
    DeleteValue delete_value;
    auto state = state_.load();
    for (size_t i = 0; i < state->buckets; ++i) {
      auto &kv = state->table[i];
      if (kv.first != Key()) {
        delete_value(kv.second);
      }
    }
    while (state) {
      auto old = state->old;
      free(state);
      state = old;
    }
  }
  Value Find(const Key &k) {
    auto state = state_.load(std::memory_order_relaxed);
    size_t b = Hash()(k) & (state->buckets - 1);
    size_t idx = 0;
    while (1) {
      auto &o = state->table[b];
      if (o.first == k) {
        return o.second;
      }
      if (o.first == Key()) break;
      b = (b + ++idx) & (state->buckets - 1);
    }
    return Value();
  }
  void Resize(int dir) {
    assert(dir == 1);
    auto old = state_.load(std::memory_order_acquire);
    size_t buckets = old->buckets * 2;
    State *state = (State*)malloc(sizeof(State) + sizeof(old->table[0]) * buckets);
    state->old = old;
    state->size = old->size;
    state->buckets = buckets;
    for (size_t i = 0; i < buckets; ++i) {
      state->table[i].first = Key();
    }
    for (size_t i = 0; i < old->buckets; ++i) {
      auto &o = old->table[i];
      if (o.first == Key()) continue;
      size_t b = Hash()(o.first) & (state->buckets - 1);
      size_t idx = 0;
      while (state->table[b].first != Key()) {
        b = (b + ++idx) & (state->buckets - 1);
      }
      state->table[b] = o;
    }
    state_.store(state, std::memory_order_release);
  }
  std::pair Add(const Key &k, const Value &value) {
    std::lock_guard lk(add_mutex_);
    auto state = state_.load(std::memory_order_acquire);
    if (state->size * 5 >= 4 * state->buckets) {
      Resize(1);
      state = state_.load(std::memory_order_relaxed);
    }
    size_t b = Hash()(k) & (state->buckets - 1);
    size_t idx = 0;
    while (1) {
      auto &o = state->table[b];
      if (o.first == k) {
        return std::make_pair(o.second, false);
      }
      if (o.first == Key()) {
        o.second = value;
        o.first = k;
        ++state->size;
        break;
      }
      b = (b + ++idx) & (state->buckets - 1);
    }
    return std::make_pair(value, true);
  }
};
Using this as the store to replace google::dense_hash_map, gives:
table: 3.96917 nano-seconds/call
table: 2.51942e+08 #/second
The Table implementation achieve the performance in cost of memory, about half bytes are wasted, and assuming that the value is never changed, deleted. And it also assumes Key() is the empty key.

TODO: open multiple method, with ambiguous resolutions.


reorder to improve cache locality

In large machine learning problems, such as parameter server, we have to calculate over samples. And for every sample, we have to fetch some variables from remote server. To improve performance, we can cache the recent used variables in local memory, and only fetch ones not cached. The problem is that, can we reorder the samples, to minimize the fetching times, and number of fetching variables?

Another case is to reorder a sparse matrix to improve the performance of matrix and vector multiplication, esp. if the sparse matrix is unchanged(or, only values are changed, the sparse structure is kept) and used many times. 

The precise problem corresponding to the best cache algorithm. I think is very hard to solve.

A greedy solution is choose the sample from not calculated samples, whose number of fetching variables(not in cache) is smallest.

In the matrix-vector example, if the sparse matrix is not randomly generated, such as user-book rating, there is a lots of people rating almost the same books, then, I guess, the performance improvement may be considerable.

In the parameter server case, the greedy algorithm works for one machine. But we still needs a algorithm to distribute the samples over machines.

代码可读性

不少程序员信奉代码可读性教条, 认为, 代码会被很多人读, 只是偶尔执行.

我的观念恰好相反. 代码会被执行很多遍, 不然, 就不必要写. 被执行很多遍的代码, 需要优化, 极尽能事地优化. 代码量大了时, 要将代码分块:

块间接口要清晰, 块内实现要高效.

被优化的代码, 可能你自己也将看不懂. 那没关系, 你只要能弄明白每一块代码的意图, 而不需要时时让每个人清楚这块代码每行的原理. 不少程序员只是机器, 会解码每行代码, 却不懂合起来的意思, 还责怪别人的代码难懂. 给定代码意图, 只写让所有程序员都理解的代码, 无谓地低效.

如果某人不清楚某块代码的意图, 那是注释的事, 你需要加注释, 以防将来你自己也忘了. 如果某人不清楚某块代码的实现原理, 那是他自己的事. 机器也不懂, 它不需要懂, 它只需要执行. 不懂某块代码的实现原理, 但知道意图, 就很容易重写, 维护问题也就不存在.

自志, 请勿对号入座.

附四年前的代码一段. 估计大部分人不会(立即)清楚这个实现原理, 但清楚意图(仅名字 Single value Linear Regression Analysis 就够了). 当然, 这段还不够高效.


class SLRA:
    '''
    SLRA: Single value Linear Regression Analysis
    usage:
      # Initial a SLRA:
      slra = SLRA(dim)
      # For every incoming vector x of dimension dim, corresponding
      # scalar result y, you can add it to the system by:
      slra.add(x, y)
      # To get the coefficients:
      slra.get()
    '''

    def __init__(self, dim, bigo=1.0e12, eps=1.0e-12):
        self.dim = dim
        self.P = [[i == j and bigo or 0.0 for i in range(dim)] for j in range(dim)]
        self.theta = [eps for i in range(dim)]

    def add(self, x, y):
        px = [sum([self.P[i][j] * x[j] for j in range(self.dim)]) for i in range(self.dim)]
        xpx1 = 1.0 + sum([x[i] * px[i] for i in range(self.dim)])
        K = [px[i] / xpx1 for i in range(self.dim)]
        delta = y - sum([x[i] * self.theta[i] for i in range(self.dim)])
        for i in range(self.dim):
            self.theta[i] += K[i] * delta
            for j in range(self.dim):
                self.P[i][j] -= K[i]*K[j]* xpx1

    def get(self):
        return self.theta